Thank you for your answer.

Well, I have already known that relation, but I was motivated by the fomula.

$$\exp\left(\sum_k a_k\cos(k\theta+\phi_k)\right)=\sum_{n=0}^\infty\frac{1}{n!}\left(\sum_{k=0}^\infty a_k\cos(k\theta+\phi_k)\right)^n$$

Although this is difficult to calculate integral as it is,
after some deformations (which possibly contain some mistakes), this will be:

$$=\sum_{n=0}^\infty\frac{1}{2^{n+1}n!}\sum_{p_0=0}^\infty\sum_{p_1=0}^\infty\cdots\sum_{p_n=0}^\infty\left(\prod_{l=0}^n a_{p_l}\right)\sum_{q_0=\pm1}\sum_{q_1=\pm1}\cdots\sum_{q_n=\pm1}e^{i\sum_{m=0}^nq_m\phi_{p_m}}\cdot e^{i\sum_{m=0}^nq_m p_m\theta}$$

This (some complicated & formal) notation can be integrated,

$$c_j:=\int_{-\pi}^{\pi}d\theta\exp\left(\sum_k a_k\cos(k\theta+\phi_k)\right)\cos j\theta$$
$$=\sum_{n=0}^\infty\frac{1}{2^{n+1}n!}\sum_{p_0=0}^\infty\sum_{p_1=0}^\infty\cdots\sum_{p_n=0}^\infty\left(\prod_{l=0}^n a_{p_l}\right)\sum_{q_0=\pm1}\sum_{q_1=\pm1}\cdots\sum_{q_n=\pm1}e^{i\sum_{m=0}^nq_m\phi_{p_m}}\cdot(\delta_{\sum_{m=0}^n q_m p_m+j}+\delta_{\sum_{m=0}^n q_m p_m-j})$$
where $\delta_r$ represents Kronecker delta.

It is difficult to make this more simple for me, but the interior of the summation seems to decrease rapidly with "well-behaved" $a_k$.

So this Fourier integral may be dominated by the smaller $n$ terms and I have managed to write down $n=0$ and $n=1$.

Then the approximated Fourier cosine series for small $k$ :
$$c_0\sim\frac{1}{2}a_0+\frac{1}{4}\left(a_0^2+\sum_{p=0}^\infty a_p^2\right)$$

$$c_1\sim\frac{1}{2}a_1\cos\phi_1+\frac{1}{2}\left(a_0a_1\cos(\phi_1)+\sum_{p=0}^{\infty}a_p a_{p+1}\cos(\phi_p-\phi_{p+1})\right)$$

As I am originally interested in lower-order terms, possibly this caluclation will be sufficient (not tested yet).