傅里叶变换（FT）： F ( ω ) = F [ f ( t ) ] = ∫ ? ∞ ∞ f ( t ) e ? j ω t d t F(\omega)=F[f(t)]=\int_{-\infty}^{ \infty}f(t)e^{-j\omega t}dt F(ω)=F[f(t)]=∫?∞ ∞f(t)e?jωtdt 傅里叶离散变换（DFT）： X ( k ) = ∑ 0 N ? 1 x ( n ) W N k n ( k = 0 , 1 , 2 , 3 ? N ? 1 ) X(k)=\sum_0^{N-1}x(n)W_N^{kn} \quad (k=0,1,2,3\cdots N-1) X(k)=0∑N−1​x(n)WNkn​(k=0,1,2,3⋯N−1) 其中 W N k n = e − j 2 π N k n W_N^{kn}=e^{-j\frac{2\pi}{N}kn} WNkn​=e−jN2π​kn，称为旋转因子。

\quad FFT是基于DFT的一种算法，将长序列的DFT分解为短序列的DFT，利用旋转因子的周期性、对称性、可约性，减少了重复运算。

对于一个多项式 A ( x ) = ∑ 0 n − 1 a i x i = a 0 + a 1 x + a 2 x 2 + ⋯ + a n − 1 x n − 1 A(x)=\sum_0^{n-1}a_ix^i=a_0+a_1x+a_2x^2+\cdots +a_{n-1}x^{n-1} A(x)=0∑n−1​ai​xi=a0​+a1​x+a2​x2+⋯+an−1​xn−1 按照下标的奇偶性把A(x)分成两部分 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 ) A(x)=(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+a_{n-1}x^{n-1}) \\ \quad\quad\; =(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+a_5x^4+a_{n-1}x^{n-2}) A(x)=(a0​+a2​x2+a4​x4+⋯+an−2​xn−2)+(a1​x+a3​x3+a5​x5+an−1​xn−1)=(a0​+a2​x2+a4​x4+⋯+an−2​xn−2)+x(a1​+a3​x2+a5​x4+an−1​xn−2) 发现奇偶部分结构相同，系数存在差异。 对于DFT来说： X ( k ) = ∑ n 为 偶 数 x ( n ) W N k n + ∑ n 为 奇 数 x ( n ) W N k n = ∑ l = 0 N / 2 − 1 x ( 2 l ) W N k 2 l + ∑ l = 0 N / 2 − 1 x ( 2 l + 1 ) W N k ( 2 l + 1 ) 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)} X(k)=n为偶数∑​x(n)WNkn​+n为奇数∑​x(n)WNkn​=l=0∑N/2−1​x(2l)WNk2l​+l=0∑N/2−1​x(2l+1)WNk(2l+1)​ 令 x 1 ( l ) = x ( 2 l ) , x 2 ( l ) = x ( 2 l + 1 ) x_1(l)=x(2l),x_2(l)=x(2l+1) x1​(l)=x(2l),x2​(l)=x(2l+1)，注意根据可约性有 W N 2 k l = W N / 2 k l W_N^{2kl}=W_{N/2}^{kl} WN2kl​=WN/2kl​ 那么 X ( k ) = ∑ l = 0 N / 2 − 1 x 1 ( l ) W N / 2 k l + W N k ∑ l = 0 N / 2 − 1 x 2 ( l ) W N / 2 k l 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} X(k)=l=0∑N/2−1​x1​(l)WN/2kl​+WNk​l=0∑N/2−1​x2​(l)WN/2kl​ 便于表示 X 1 ( k ) = ∑ l = 0 N / 2 − 1 x 1 ( l ) W N / 2 k l X 2 ( k ) = ∑ l = 0 N / 2 − 1 x 2 ( l ) W N / 2 k l 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} X1​(k)=l=0∑N/2−1​x1​(l)WN/2kl​X2​(k)=l=0∑N/2−1​x2​(l)WN/2kl​ 均以N/2为周期 利用旋转因子对称性 W N m + N / 2 = − W N m \quad W_N^{m+N/2}=-W_N^{m} WNm+N/2​=−WNm​

当 k = 0 , 1 , 2 , ⋯   , N 2 − 1 k=0,1,2,\cdots,\frac{N}{2}-1 k=0,1,2,⋯,2N​−1时: X ( k ) = X 1 ( k ) + W N k X 2 ( k ) X ( k + N 2 ) = X 1 ( k + N 2 ) + W N k + N / 2 X 2 ( k + N 2 ) = X 1 ( k ) − W N k X 2 ( k ) X(k)=X_1(k)+W_N^kX_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^kX_2(k) X(k)=X1​(k)+WNk​X2​(k)X(k+2N​)=X1​(k+2N​)+WNk+N/2​X2​(k+2N​)=X1​(k)−WNk​X2​(k) 可以看到，只需要一半的数据即可获得整个序列的离散傅里叶变换。 基础的，可以通过递归分治的方法计算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按照奇偶区分的数组。 计算好的对应傅里叶变换值直接存到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;
			}
		}
	}
}