從頻域計算信號回圈譜

回圈譜的計算最大的特點就是復雜，在知網上，嵌入式系統上的快速回圈譜計算的文章也非常多，不過在理論分析階段，只要能加快仿真計算就夠了，

標準的直接從時域信號計算回圈自相關的公式為：
R x α ( τ ) = lim ? Δ t → ∞ 2 Δ t ∫ ? Δ t / 2 Δ t / 2 x ( t + τ / 2 ) x ? ( t ? τ / 2 ) e ? j 2 π α t d t R_{x}^{\alpha}(\tau)=\lim_{\Delta{t}\rightarrow\infty}\frac{2}{\Delta{t}} \int_{-\Delta{t}/2}^{\Delta{t}/2}x(t+\tau/2)x^*(t-\tau/2)e^{-j2\pi\alpha t}dt Rxα?(τ)=Δt→∞lim?Δt2?∫?Δt/2Δt/2?x(t+τ/2)x?(t?τ/2)e?j2παtdt
其中 α \alpha α為回圈頻率，我們在計算機處理中，自然都要離散化，將 t t t與 n n n對應， τ \tau τ與 m m m對應，上式就轉變為
R x α ( m ) = 1 N ∑ n = 0 N ? 1 x ( n + m ) x ? ( n ) exp ? ( ? j 2 π α n ) exp ? ( ? j π α m ) R_x^{\alpha}(m)=\frac{1}{N}\sum_{n=0}^{N-1}x(n+m)x^*(n)\exp(-j2\pi\alpha{n})\exp(-j\pi\alpha{m}) Rxα?(m)=N1?n=0∑N?1?x(n+m)x?(n)exp(?j2παn)exp(?jπαm)
而回圈譜就是回圈自相關的傅里葉變換，離散時可以用DFT：
S x α ( k ) = DFT [ R x α ( m ) ] S_x^{\alpha}(k)=\text{DFT}\left[R_x^{\alpha}(m)\right] Sxα?(k)=DFT[Rxα?(m)]

這里最大的問題是要對信號求相關，浪費大量的時間，但時域相卷對應頻域相乘，時域的相關也可以用DFT變成頻域的逐點相乘，我們對第2式做DFT，可得：
S x α ( k ) = 1 N ∑ m = 0 N ? 1 ∑ n = 0 N ? 1 x ( n + m ) x ? ( n ) exp ? ( ? j 2 π α n ) exp ? ( ? j π α m ) exp ? ( ? j 2 π k m N ) S_x^{\alpha}(k)=\frac{1}{N}\sum_{m=0}^{N-1} \sum_{n=0}^{N-1}x(n+m)x^*(n)\exp(-j2\pi\alpha{n})\exp(-j\pi\alpha{m}) \exp(-j\frac{2\pi km}{N}) Sxα?(k)=N1?m=0∑N?1?n=0∑N?1?x(n+m)x?(n)exp(?j2παn)exp(?jπαm)exp(?jN2πkm?)
整理拆分得到：
S x α ( k ) = 1 N [ ∑ m = 0 N ? 1 x ( n + m ) e ? j π α ( n + m ) exp ? ( ? j 2 π k ( m + n ) N ) ] [ ∑ n = 0 N ? 1 x ( n ) e j π α n exp ? ( ? j 2 π k n 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]^* Sxα?(k)=N1?[m=0∑N?1?x(n+m)e?jπα(n+m)exp(?jN2πk(m+n)?)][n=0∑N?1?x(n)ejπαnexp(?jN2πkn?)]?
只要 N α / 2 N\alpha/2 Nα/2為正數，就有
S x α ( k ) = 1 N X ( k + N α / 2 ) X ? ( k ? N α / 2 ) S_x^{\alpha}(k)=\frac{1}{N}X(k+N\alpha/2)X^*(k-N\alpha/2) Sxα?(k)=N1?X(k+Nα/2)X?(k?Nα/2)
這樣就可以很方便地直接從頻域獲得回圈譜，相比自相關 O ( n 2 ) O(n^2) O(n2)的計算量，做FFT需要 O ( n log ? ( n ) ) O\left(n\log(n)\right) O(nlog(n))，頻域相乘只要 O ( n ) O(n) O(n)，

代碼掛在Github上：CyclicSpectrumPlot

在實作時，我在做FFT時多補一倍的零，使得 α / 2 \alpha/2 α/2可以更簡便離散化計算，對應的，我的程式會引入更大的柵瓣，