FFT的源码含义

Posted 2023-02-26

参考技术A

在C环境下的源码
源码（1）：
注：亲测，这个版本无法运行，作者删改了重要内容 请参考源码（2） //快速傅立叶变换// 入口参数：// l: l=0, 傅立叶变换;l=1, 逆傅立叶变换// il: il=0,不计算傅立叶变换或逆变换模和幅角；il=1,计算模和幅角// n: 输入的点数，为偶数，一般为32，64，128，...,1024等// k: 满足n=2^k(k>0),实质上k是n个采样数据可以分解为偶次幂和奇次幂的次数// pr[]: l=0时，存放N点采样数据的实部// l=1时, 存放傅立叶变换的N个实部// pi[]: l=0时，存放N点采样数据的虚部// l=1时, 存放傅立叶变换的N个虚部//// 出口参数：// fr[]: l=0, 返回傅立叶变换的实部// l=1, 返回逆傅立叶变换的实部// fi[]: l=0, 返回傅立叶变换的虚部// l=1, 返回逆傅立叶变换的虚部// pr[]: il=1,i=0 时，返回傅立叶变换的模// il=1,i=1 时，返回逆傅立叶变换的模// pi[]: il=1,i=0 时，返回傅立叶变换的辐角// il=1,i=1 时，返回逆傅立叶变换的辐角void fft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il)int it,m,is,i,j,nv,l0;double p,q,s,vr,vi,poddr,poddi;for(it=0;it<=n-1;m=it++)is=0;for(i=0;i<=k-1;i++)j=m/2;is=2*is+(m-2*j);m=j;fr[it]=pr[is];fi[it]=pi[is];//----------------------------pr[0]=1.0;pi[0]=0.0;p=6.283185306/n;pr[1]=cos(p);pi[1]=-sin(p);if (l)pi[1]=-pi[1];for(i=2;i<=n-1;i++)p=pr[i-1]*pr[1];q=pi[i-1]*pi[1];s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);pr=p-q;pi=s-p-q;for(it=0;it<=n-2;it+=2)vr=fr[it];vi=fi[it];fr[it]=vr+fr[it+1];fi[it]=vi+fi[it+1];fr[it+1]=vr-fr[it+1];fi[it+1]=vi-fi[it+1];m=n/2;nv=2;for(l0=k-2;l0>=0;l0--)m/=2;nv<<=1;for(it=0;it<=(m-1)*nv;it+=nv)for(j=0;j<=(nv/2)-1;j++)p=pr[m*j]*fr[it+j+nv/2];q=pi[m*j]*fi[it+j+nv/2];s=pr[m*j]+pi[m*j];s*=(fr[it+j+nv/2]+fi[it+j+nv/2]);poddr=p-q;poddi=s-p-q;fr[it+j+nv/2]=fr[it+j]-poddr;fi[it+j+nv/2]=fi[it+j]-poddi;fr[it+j]+=poddr;fi[it+j]+=poddi;if(l)for(i=0;i<=n-1;fr/=n,fi[i++]/=n);if(il)for(i=0;i<=n-1;i++)pr=sqrt(fr*fr+fi*fi);if(fabs(fr)<0.000001*fabs(fi))pi=fi*fr>0?90.0-90.0;elsepi=atan(fi/fr)*360.0/6.283185306;return;源码（2）
ps:可以运行的 // The following line must be defined before including math.h to correctly define M_PI#define _USE_MATH_DEFINES#include <math.h>#include <stdio.h>#include <stdlib.h>#define PI M_PI /* pi to machine precision, defined in math.h */#define TWOPI (2.0*PI)/*FFT/IFFT routine. (see pages 507-508 of Numerical Recipes in C)Inputs:data[] : array of complex* data points of size 2*NFFT+1.data[0] is unused,* the n'th complex number x(n), for 0 <= n <= length(x)-1, is stored as:data[2*n+1] = real(x(n))data[2*n+2] = imag(x(n))if length(Nx) < NFFT, the remainder of the array must be padded with zerosnn : FFT order NFFT. This MUST be a power of 2 and >= length(x).isign: if set to 1,computes the forward FFTif set to -1,computes Inverse FFT - in this case the output values haveto be manually normalized by multiplying with 1/NFFT.Outputs:data[] : The FFT or IFFT results are stored in data, overwriting the input.*/void four1(double data[], int nn, int isign)int n, mmax, m, j, istep, i;double wtemp, wr, wpr, wpi, wi, theta;double tempr, tempi;n = nn << 1;j = 1;for (i = 1; i < n; i += 2) if (j > i) tempr = data[j]; data[j] = data[i]; data[i] = tempr;tempr = data[j+1]; data[j+1] = data[i+1]; data[i+1] = tempr;m = n >> 1;while (m >= 2 && j > m) j -= m;m >>= 1;j += m;mmax = 2;while (n > mmax) istep = 2*mmax;theta = TWOPI/(isign*mmax);wtemp = sin(0.5*theta);wpr = -2.0*wtemp*wtemp;wpi = sin(theta);wr = 1.0;wi = 0.0;for (m = 1; m < mmax; m += 2) for (i = m; i <= n; i += istep) j =i + mmax;tempr = wr*data[j] - wi*data[j+1];tempi = wr*data[j+1] + wi*data[j];data[j] = data[i] - tempr;data[j+1] = data[i+1] - tempi;data[i] += tempr;data[i+1] += tempi;wr = (wtemp = wr)*wpr - wi*wpi + wr;wi = wi*wpr + wtemp*wpi + wi;mmax = istep;在C++环境下的源码
bool FFT(complex<double> * TD, complex<double> * FD, int r)//一维快速Fourier变换。//complex<double> * TD ——指向时域数组的指针; complex<double> * FD ——指向频域数组的指针; r ——2的幂数，即迭代次数LONG count; // Fourier变换点数int i,j,k; // 循环变量int bfsize,p; // 中间变量double angle; // 角度complex<double> *W,*X1,*X2,*X;count = 1 << r; // 计算Fourier变换点数为1左移r位W = new complex<double>[count / 2];X1 = new complex<double>[count];X2 = new complex<double>[count]; // 分配运算所需存储器// 计算加权系数（旋转因子w的i次幂表）for(i = 0; i < count / 2; i++)angle = -i * PI * 2 / count;W[ i ] = complex<double> (cos(angle), sin(angle));// 将时域点写入X1memcpy(X1, TD, sizeof(complex<double>) * count);// 采用蝶形算法进行快速Fourier变换for(k = 0; k < r; k++)for(j = 0; j < 1 << k; j++)bfsize = 1 << (r-k);for(i = 0; i < bfsize / 2; i++)p = j * bfsize;X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2] * W[i * (1<<k)];X2[i + p + bfsize / 2] = X1[i + p] - X1[i + p + bfsize / 2] * W[i * (1<<k)];X = X1;X1 = X2;X2 = X;// 重新排序for(j = 0; j < count; j++)p = 0;for(i = 0; i < r; i++)if (j&(1<<i))p+=1<<(r-i-1);FD[j]=X1[p];// 释放内存delete W;delete X1;delete X2;return true;在Matlab环境下的源码
function X=myfft(x)
N=length(x);
if N==1
X=x;
return;
end
%myfft函数 用递归实现
t=log2(N);
t1=floor(t);
t2=ceil(t);
if t1~=t2; %若x的长度N不为2的整数次幂，则补0至最接近的2的整数次幂
x=[x zeros(1,2^t2-N)];
N=2^t2;
end
w0=exp(-1i*2*pi/N);
X=zeros(1,N);
n=1:N/2;
xo(n)=x(2*n-1);
xe(n)=x(2*n);
XO=myfft(xo); %递归调用
XE=myfft(xe);
for n=1:N/2
X(n)=XO(n)+XE(n)*(w0^(n-1));
X(n+N/2)=XO(n)-XE(n)*(w0^(n-1));
end