C++实现快速傅里叶变换_c++傅里叶变换程序

基础公式

傅里叶变换（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\cdots N-1)$$
其中$W_N^{kn}=e^{-j\frac{2\pi }{N}kn}$，称为旋转因子。

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

FFT原理

对于一个多项式
$$A(x)=\sum _0^{n-1}{a}_i{x}^i=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}$$
按照下标的奇偶性把A(x)分成两部分
$$\begin{gathered} A(x)=(a_0+a_2{x}^2+a_4{x}^4+\cdots +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+\cdots +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,\cdots ,\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;
};
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

递归法

定义两个复数类数组，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;
		}
	}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

验证代码，设置采样点数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();
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

结果：

可以看到频率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;
			}
		}
	}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27