C语言、Matlab实现FFT几种编程实例..

C语言、MATLAB实现FFT几种方法

总结前人经验，仅供参考

///一、

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////c语言程序////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// #include

#include

#include

#define PI 3.1415926535897932384626433832795028841971 //定义圆周率值#define FFT_N 128 //定义福利叶变换的点数

struct compx {float real,imag;}; //定义一个复数结构struct compx s[FFT_N]; //FFT输入和输出：从S[1]开始存放，根据大小自己定义/******************************************************************* 函数原型：struct compx EE(struct compx b1,struct compx b2)

函数功能：对两个复数进行乘法运算

输入参数：两个以联合体定义的复数a,b

输出参数：a和b的乘积，以联合体的形式输出

*******************************************************************/ struct compx EE(struct compx a,struct compx b)

{

struct compx c;

c.real=a.real*b.real-a.imag*b.imag;

c.imag=a.real*b.imag+a.imag*b.real;

return(c);

}

/*****************************************************************

函数原型：void FFT(struct compx *xin,int N)

函数功能：对输入的复数组进行快速傅里叶变换（FFT）

输入参数：*xin复数结构体组的首地址指针，struct型

*****************************************************************/ void FFT(struct compx *xin)

{

int f,m,nv2,nm1,i,k,l,j=0;

struct compx u,w,t;

nv2=FFT_N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法nm1=FFT_N-1;

for(i=0;i

{

if(i

{

t=xin[j];

xin[j]=xin[i];

xin[i]=t;

}

k=nv2; //求j的下一个倒位序

while(k<=j) //如果k<=j,表示j的最高位为1

{

j=j-k; //把最高位变成0

k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0

}

j=j+k; //把0改为1

}

{

int le,lei,ip; //FFT运算核，使用蝶形运算完成FFT运算f=FFT_N;

for(l=1;(f=f/2)!=1;l++) //计算l的值，即计算蝶形级数;

for(m=1;m<=l;m++) // 控制蝶形结级数

{ //m表示第m级蝶形，l为蝶形级总数l=log（2）N le=2<<(m-1); //le蝶形结距离，即第m级蝶形的蝶形结相距le点lei=le/2; //同一蝶形结中参加运算的两点的距离

u.real=1.0; //u为蝶形结运算系数，初始值为1

u.imag=0.0;

w.real=cos(PI/lei); //w为系数商，即当前系数与前一个系数的商w.imag=-sin(PI/lei);

for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结{

for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算，即计算系数相同蝶形结

{

ip=i+lei; //i，ip分别表示参加蝶形运算的两个节点

t=EE(xin[ip],u); //蝶形运算，详见公式

xin[ip].real=xin[i].real-t.real;

xin[ip].imag=xin[i].imag-t.imag;

xin[i].real=xin[i].real+t.real;

xin[i].imag=xin[i].imag+t.imag;

}

u=EE(u,w); //改变系数，进行下一个蝶形运算

}

}

}

}

/************************************************************

函数原型：void main()

函数功能：测试FFT变换，演示函数使用方法

输入参数：无

输出参数：无

************************************************************/

void main()

{

int i;

for(i=0;i

{

s[i].real=sin(2*3.141592653589793*i/FFT_N); //实部为正弦波FFT_N点采样，赋值为1

s[i].imag=0; //虚部为0

}

FFT(s); //进行快速福利叶变换

for(i=0;i

while(1);

}

%////二、

%/////////////////////////////////////////////////////////////////////////////////////////////////////////////

%///////////////////////////////////////////////////////////////////////////////////////////////////////////

%////////////////////////////////MATLAB仿真信号源的源程序::

Clear;

Clc;

t=O:O.01:3;

yl=100*sin(pi/3*t);

n=l;

for t=-O:O.01:10;

y2(1,n)=-61.1614*exp(-0.9*t);

n=n+;

end

min(y2)

y=[yl，y2];

figure(1);

plot(y); %funboxing(O.001+1/3)

%////////////////////////

%/////////快速傅里叶变换matlab程序:

%////////////////////////clc;

clear;

clf;

N=input('Node number')

T=input('cai yang jian ge')

f=input('frenquency')

choise=input('add zero or not? 1/0 ')

n=0:T:(N-1)*T; %采样点

k=0:N-1;

x=sin(2*pi*f*n);

if choise==1

e=input('Input the number of zeros!')

x=[x zeros(1,e)]

N=N+e;

else

end %加零

k=0:N-1; %给k重新赋值，因为有可能出现加零状况bianzhi=bi2de(fliplr(de2bi(k,length(de2bi(N-1)))))+1;%利用库函数进行变址运算for l=1:N

X(l)=x(bianzhi(l));%将采样后的值按照变址运算后的顺序放入X矩阵中end

d1=1;

for m=1:log2(N)

d2=d1; %做蝶形运算的两个数之间的距离

d1=d1*2; %同一级之下蝶形结之间的距离

W=1; %蝶形运算系数的初始值

dw=exp(-j*pi/d2); %蝶形运算系数的增加量

for t=1:d2 %

for i=t:d1:N

i1=i+d2;

if(i1>N)break; %判断是否超出范围

else

p=X(i1)*W;

X(i1)=X(i)-p;

X(i)=X(i)+p; %蝶形运算

end

end

W=W*dw; %蝶形运算系数的变化endend

subplot(2,2,1);

t=0:0.0000001:N*T;

plot(t,sin(2*pi*f*t)); %画原曲线

subplot(2,2,2);

stem(k,x); %画采样后的离散信号

subplot(2,2,3);

stem(k,abs(X)/max(abs(X))); %画自己的fft的运算结果

subplot(2,2,4);

stem(k,abs(fft(x))/max(abs(fft(x)))); %调用系统的函数绘制结果，与自己的结果进行比较

//////三、

///////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////FFT的C语言算法实现

////////////程序如下：

/************FFT***********/

#include

#include

#include

#define N 1000

typedef struct

{

double real;

double img;

}complex;

void fft(); /*快速傅里叶变换*/

void ifft(); /*快速傅里叶逆变换*/

void initW();

void change();

void add(complex ,complex ,complex *); /*复数加法*/

void mul(complex ,complex ,complex *); /*复数乘法*/

void sub(complex ,complex ,complex *); /*复数减法*/

void divi(complex ,complex ,complex *);/*复数除法*/ void output(); /*输出结果*/

complex x[N], *W;/*输出序列的值*/

int size_x=0;/*输入序列的长度，只限2的N次方*/

double PI;

int main()

{

int i,method;

system("cls");

PI=atan(1)*4;

printf("Please input the size of x:\n");

/*输入序列的长度*/

scanf("%d",&size_x);

printf("Please input the data in x[N]:(such as:5 6)\n"); /*输入序列对应的值*/

for(i=0;i

scanf("%lf %lf",&x[i].real,&x[i].img);

initW();

/*选择FFT或逆FFT运算*/

printf("Use FFT(0) or IFFT(1)?\n");

scanf("%d",&method);

if(method==0)

fft();

else

ifft();

output();

return 0;

}

/*进行基-2 FFT运算*/

void fft()

{

int i=0,j=0,k=0,l=0;

complex up,down,product;

change();

for(i=0;i< log(size_x)/log(2) ;i++)

{

l=1<

for(j=0;j

{

for(k=0;k

{

mul(x[j+k+l],W[size_x*k/2/l],&product);

add(x[j+k],product,&up);

sub(x[j+k],product,&down);

x[j+k]=up;

x[j+k+l]=down;

}

}

}

}

void ifft()

{

int i=0,j=0,k=0,l=size_x;

complex up,down;

for(i=0;i< (int)( log(size_x)/log(2) );i++) /*蝶形运算*/ {

l/=2;

for(j=0;j

{

for(k=0;k

{

add(x[j+k],x[j+k+l],&up);

up.real/=2;up.img/=2;

sub(x[j+k],x[j+k+l],&down);

down.real/=2;down.img/=2;

divi(down,W[size_x*k/2/l],&down);

x[j+k]=up;

x[j+k+l]=down;

}

}

}

change();

}

void initW()

{

int i;

W=(complex *)malloc(sizeof(complex) * size_x); for(i=0;i

{

W[i].real=cos(2*PI/size_x*i);

W[i].img=-1*sin(2*PI/size_x*i);

}

}

void change()

{

complex temp;

unsigned short i=0,j=0,k=0;

double t;

for(i=0;i

{

k=i;j=0;

t=(log(size_x)/log(2));

while( (t--)>0 )

{

j=j<<1;

j|=(k & 1);

k=k>>1;

}

if(j>i)

{

temp=x[i];

x[i]=x[j];

x[j]=temp;

}

}

}

void output() /*输出结果*/

{

int i;

printf("The result are as follows\n");

for(i=0;i

{

printf("%.4f",x[i].real);

if(x[i].img>=0.0001)

printf("+%.4fj\n",x[i].img);

else if(fabs(x[i].img)<0.0001)

printf("\n");

else

printf("%.4fj\n",x[i].img);

}

}

void add(complex a,complex b,complex *c) {

c->real=a.real+b.real;

c->img=a.img+b.img;

}

void mul(complex a,complex b,complex *c) {

c->real=a.real*b.real - a.img*b.img;

c->img=a.real*b.img + a.img*b.real;

}

void sub(complex a,complex b,complex *c) {

c->real=a.real-b.real;

c->img=a.img-b.img;

}

void divi(complex a,complex b,complex *c)

{

c->real=( a.real*b.real+a.img*b.img )/(

b.real*b.real+b.img*b.img);

c->img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);

}

/////四、

///////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// %//////选带傅里叶变换(zoom-fft) MATLAB

%移频（将选带的中心频率移动到零频）

%数字低通滤波器（防止频率混叠）

%重新采样（将采样的数据再次间隔采样，间隔的数据取决于分析的带宽，就是放大倍数）

%复FFT （由于经过了移频，所以数据不是实数了）

%频率调整（将负半轴的频率成分移到正半轴）

function [f, y] = zfft(x, fi, fa, fs)

% x为采集的数据

% fi为分析的起始频率

% fa为分析的截止频率

% fs为采集数据的采样频率

% f为输出的频率序列

% y为输出的幅值序列(实数）

f0 = (fi + fa) / 2; %中心频率

N = length(x); %数据长度

r = 0:N-1;

b = 2*pi*f0.*r ./ fs;

x1 = x .* exp(-1j .* b); %移频

bw = fa - fi;

B = fir1(32, bw / fs); %滤波截止频率为0.5bw

x2 = filter(B, 1, x1);

c = x2(1:floor(fs/bw):N); %重新采样

N1 = length(c);

f = linspace(fi, fa, N1);

y = abs(fft(c)) ./ N1 * 2;

y = circshift(y, [0, floor(N1/2)]); %将负半轴的幅值移过来

end

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////

%上边四程序测试应用实例：

fs = 2048;

T = 100;

t = 0:1/fs:T;

x = 30 * cos(2*pi*110.*t) + 30 * cos(2*pi*111.45.*t) + 25*cos(2*pi*112.3*t) + 48*c os(2*pi*113.8.*t)+50*cos(2*pi*114.5.*t);

[f, y] = zfft(x, 109, 115, fs);

plot(f, y);

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

图片效果：