FFT源代码



源码表示编辑

在C环境下的源码

源码（1）：

注：亲测，这个版本无法运行，作者删改了重要内容^[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;

else

pi=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 zeros

nn : FFT order NFFT. This MUST be a power of 2 and >= length(x).

isign: if set to 1,

computes the forward FFT

if set to -1,

computes Inverse FFT - in this case the output values have

to 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));

}

// 将时域点写入X1

memcpy(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)

%myfft函数 用递归实现

N=length(x);

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(-j*2*pi/N);

X=zeros(1,N);

if N==2

X(1)=x(1)+x(2);

X(2)=x(1)-x(2);

else

n=1:N/2;

xe(n)=x(2*n-1);

xo(n)=x(2*n);

XE=myfft(xe); %递归调用

XO=myfft(xo);

for n=1:N/2

X(n)=XE(n)+XO(n)*(w0^(n-1));

X(n+N/2)=XE(n)-XO(n)*(w0^(n-1));

end

end