基礎公式

傅里葉變換（FT）：
$$F(\omega )=F[f(t)]={\int }_{?\infty }^{+\infty }f(t)e^{?j\omega t}dt$$
離散傅里葉變換（DFT）：
$$X(k)=\sum _0^{N?1}x(n)W_N^{kn}(k=0,1,2,3?N?1)$$
其中$W_N^{kn}=e^{?j\frac{2\pi }{N}kn}$，稱為旋轉因子。

\quadFFT是基于DFT的一種算法，將長序列的DFT分解為短序列的DFT，利用旋轉因子的周期性、對稱性、可約性，減少了重復運算。

FFT原理

對于一個多項式
$$A(x)=\sum _0^{n?1}{a}_i{x}^i=a_0+a_{1}x+a_2{x}^2+?+a_{n?1}{x}^{n?1}$$
按照下標的奇偶性把A(x)分成兩部分
$$\begin{gathered} A(x)=(a_0+a_2{x}^2+a_4{x}^4+?+a_{n?2}{x}^{n?2})+(a_{1}x+a_3{x}^3+a_5{x}^5+a_{n?1}{x}^{n?1}) \\ =(a_0+a_2{x}^2+a_4{x}^4+?+a_{n?2}{x}^{n?2})+x(a_1+a_3{x}^2+a_5{x}^4+a_{n?1}{x}^{n?2}) \end{gathered}$$
發現奇偶部分結構相同，系數存在差異。
對于DFT來說：
$$\begin{gathered} X(k)=\sum _{n為偶數}x(n)W_N^{kn}+\sum _{n為奇數}x(n)W_N^{kn} \\ =\sum _{l=0}^{N/2?1}x(2l)W_N^{k2l}+\sum _{l=0}^{N/2?1}x(2l+1)W_N^{k(2l+1)} \end{gathered}$$
令$x_1(l)=x(2l),x_2(l)=x(2l+1)$，注意根據可約性有$W_N^{2kl}=W_{N/2}^{kl}$
那么$$X(k)=\sum _{l=0}^{N/2?1}{x}_1(l)W_{N/2}^{kl}+W_N^k\sum _{l=0}^{N/2?1}{x}_2(l)W_{N/2}^{kl}$$
便于表示$$\begin{gathered} X_1(k)=\sum _{l=0}^{N/2?1}{x}_1(l)W_{N/2}^{kl} \\ {X}_2(k)=\sum _{l=0}^{N/2?1}{x}_2(l)W_{N/2}^{kl} \end{gathered}$$
均以N/2為周期
利用旋轉因子對稱性$W_N^{m+N/2}=?W_N^m$

當$k=0,1,2,?,\frac{N}{2}?1$時:
$$\begin{gathered} X(k)=X_1(k)+W_N^k{X}_2(k) \\ X(k+\frac{N}{2})=X_1(k+\frac{N}{2})+W_N^{k+N/2}{X}_2(k+\frac{N}{2})=X_1(k)?W_N^k{X}_2(k) \end{gathered}$$
可以看到，只需要一半的數據即可獲得整個序列的離散傅里葉變換。
基礎的，可以通過遞歸分治的方法計算DFT；進階的，通過遞歸發現規律，優化成蝶形算法進行計算。

下面先介紹遞歸法

復數類

因為涉及復數的運算，所以可以先設計一個復數類。
復數類應該包含復數的實部和虛部，重載運算符，使其支持加減法、乘法和求模運算。

class Complex{//復數類 public:Complex(double re=0,double im=0):real{re},image{im}{}Complex(const Complex&cp):real{cp.real},image{cp.image}{}Complex operator +(const Complex cp){return Complex(real+cp.real,image+cp.image);}Complex operator -(const Complex cp){return Complex(real-cp.real,image-cp.image);}Complex operator *(const Complex cp){double re=real*cp.real-image*cp.image;double im=image*cp.real+real*cp.image;return Complex(re,im);}double model(){return sqrt(real*real+image*image);}double real=0;double image=0; };

遞歸法

定義兩個復數類數組，resource代表數據源，result代表對resource按照奇偶區分的數組。

void Page_Tool::fftIterate(Complex *resource, Complex *result, int n) {if(n>1) { //迭代結束標志for(int i=0; i<n; i+=2) {result[i/2]=resource[i];result[i/2+n/2]=resource[i+1];}fftIterate( result,resource,n/2);fftIterate( result+n/2,resource,n/2 );for(int i=0;i<n/2;i++) {Complex omega(cos(2*PI*i/double(n)),-sin(2*PI*i/double(n)));Complex temp=omega*result[i+n/2];resource[i]=result[i]+temp;resource[i+n/2]=result[i]-temp;}} }

驗證代碼，設置采樣點數N=32768，采樣頻率Fs=32768，給定兩個不同頻率和幅值的正弦波，進行傅里葉變換，查看結果。注意，這里的分辨率為N/Fs=1Hz，波形頻率分量不滿足該分辨率時會存在頻率泄露現象。最后的結果需要進行頻譜變換才能代表實際的頻率及幅值。

void Page_Tool::on_pushButton_2_clicked() {//一些畫圖的代碼不用管ui->widget_1->m_plot->clearGraphs();ui->widget_2->m_plot->clearGraphs();ui->widget_1->addGraph();ui->widget_2->addGraph();int N=32768;//采樣點數double Fs=32768;//采樣頻率double W1=500;//波形1頻率double W2=10000;//波形2頻率Complex resource[N],result[N];for(int k=0; k<N; k++) {resource[k]=Complex(20*sin(2*PI*W1*k/Fs)+100*sin(2*PI*W2*k/Fs),0);//兩個正弦波合成信號}for(int i=0;i<N;i++){ui->widget_1->m_plot->graph(0)->addData(i/Fs,resource[i].real);}fftIterate(resource,result,N);//驗證fft代碼for(int i=0;i<N/2+1;i++){ui->widget_2->m_plot->graph(0)->addData(i*Fs/N,resource[i].model()*2/N);}ui->widget_1->m_plot->replot();ui->widget_2->m_plot->replot(); }

結果：

可以看到頻率500Hz、幅值20的正弦波和頻率10000HZ、幅值100的正弦波都被表示出來了。

蝶形法

原理大家可以搜一搜，這里給出參考代碼。效果和上面是一樣的。原始數據要組成2的整數倍，不滿足的可以補0。

void Page_Tool::fftButterFly(Complex *resource, Complex *result, int n) {int M=int(log2(n));for(int i=0;i<n;i++){int k=0;for(int j=0;j<M;j++){if(((i>>j)&1)==1){k+=1<<(M-j-1);}}result[k]=resource[i];//取二進制倒碼}for(int L=1;L<=M;L++){int B=int(pow(2,L-1));//間隔int k=int(pow(2,M-L));//增量for(int i=0;i<k;i++){//蝶形運算的次數for(int j=0;j<B;j++){//每個蝶形運算的乘法次數int index=j+2*B*i;//數組下標int p=j*k;Complex omega(cos(2*PI*p/double(n)),-sin(2*PI*p/double(n)));Complex temp=omega*result[index+B];result[index+B]=result[index]-temp;result[index]=result[index]+temp;}}} }

總結

以上是生活随笔為你收集整理的C++实现快速傅里叶变换的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。