从频域计算信号循环谱

循环谱的计算最大的特点就是复杂，在知网上，嵌入式系统上的快速循环谱计算的文章也非常多。不过在理论分析阶段，只要能加快仿真计算就够了。

标准的直接从时域信号计算循环自相关的公式为：
$$R_x^{\alpha }(\tau )=\underset{\Delta t\to \infty }{\lim }\frac{2}{\Delta t}{\int }_{-\Delta t/2}^{\Delta t/2}x(t+\tau /2)x^{\ast }(t-\tau /2)e^{-j2\pi \alpha t}dt$$
其中$\alpha$为循环频率。我们在计算机处理中，自然都要离散化，将$t$与$n$对应，$\tau$与$m$对应，上式就转变为
$$R_x^{\alpha }(m)=\frac{1}{N}\sum _{n=0}^{N-1}x(n+m)x^{\ast }(n)\exp (-j2\pi \alpha n)\exp (-j\pi \alpha m)$$
而循环谱就是循环自相关的傅里叶变换，离散时可以用DFT：
$$S_x^{\alpha }(k)=\text{DFT}\left[R_x^{\alpha }(m)\right]$$

这里最大的问题是要对信号求相关，浪费大量的时间。但时域相卷对应频域相乘，时域的相关也可以用DFT变成频域的逐点相乘。我们对第2式做DFT，可得：
$$S_x^{\alpha }(k)=\frac{1}{N}\sum _{m=0}^{N-1}\sum _{n=0}^{N-1}x(n+m)x^{\ast }(n)\exp (-j2\pi \alpha n)\exp (-j\pi \alpha m)\exp (-j\frac{2\pi km}{N})$$
整理拆分得到：
$$S_x^{\alpha }(k)=\frac{1}{N}\left[\sum _{m=0}^{N-1}x(n+m)e^{-j\pi \alpha (n+m)}\exp \left(-j\frac{2\pi k(m+n)}{N}\right)\right]{\left[\sum _{n=0}^{N-1}x(n)e^{j\pi \alpha n}\exp \left(-j\frac{2\pi kn}{N}\right)\right]}^{\ast }$$
只要$N\alpha /2$为正数，就有
$$S_x^{\alpha }(k)=\frac{1}{N}X(k+N\alpha /2)X^{\ast }(k-N\alpha /2)$$
这样就可以很方便地直接从频域获得循环谱，相比自相关$O(n^2)$的计算量，做FFT需要$O\left(n\log (n)\right)$，频域相乘只要$O(n)$。

代码挂在Github上：CyclicSpectrumPlot

在实现时，我在做FFT时多补一倍的零，使得$\alpha /2$可以更简便离散化计算。对应的，我的程序会引入更大的栅瓣。