C++实现快速傅里叶变换

阿礼先森

已于 2022-08-17 08:38:12 修改

阅读量4.3k

点赞数 3

分类专栏：控制理论文章标签： c++ 线性代数矩阵 qt matlab

于 2022-05-07 11:15:25 首次发布

本文链接：https://blog.csdn.net/jiushigedidi/article/details/124613145

版权

控制理论专栏收录该内容

1 篇文章 0 订阅

订阅专栏

基础公式

傅里叶变换（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} \quad (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_ix^i=a_0+a_1x+a_2x^2+\cdots +a_{n-1}x^{n-1}$
按照下标的奇偶性把A(x)分成两部分
$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})$
发现奇偶部分结构相同，系数存在差异。
对于DFT来说：
$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_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}$
便于表示 $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}$
均以N/2为周期
利用旋转因子对称性 $\quad W_N^{m+N/2}=-W_N^{m}$

当 $k=0,1,2,\cdots,\frac{N}{2}-1$ 时:
$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)$
可以看到，只需要一半的数据即可获得整个序列的离散傅里叶变换。
基础的，可以通过递归分治的方法计算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;
			}
		}
	}
}