完整word版C语言Matlab实现FFT几种编程实例.docx
《完整word版C语言Matlab实现FFT几种编程实例.docx》由会员分享，可在线阅读，更多相关《完整word版C语言Matlab实现FFT几种编程实例.docx（14页珍藏版）》请在冰豆网上搜索。
完整word版C语言Matlab实现FFT几种编程实例
C语言、MATLAB实现FFT几种方法
总结前人经验，仅供参考
///一、
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////c语言程序//////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
#include
#include
#include
#definePI3.1415926535897932384626433832795028841971//定义圆周率值
#defineFFT_N128//定义福利叶变换的点数
structcompx{floatreal,imag;};//定义一个复数结构
structcompxs[FFT_N];//FFT输入和输出：
从S[1]开始存放，根据大小自己定义
/*******************************************************************
函数原型：
structcompxEE（structcompxb1,structcompxb2）
函数功能：
对两个复数进行乘法运算
输入参数：
两个以联合体定义的复数a,b
输出参数：
a和b的乘积，以联合体的形式输出
*******************************************************************/
structcompxEE（structcompxa,structcompxb）
{
structcompxc;
c.real=a.real*b.real-a.imag*b.imag;
c.imag=a.real*b.imag+a.imag*b.real;
return（c）;
}
/*****************************************************************
函数原型：
voidFFT（structcompx*xin,intN）
函数功能：
对输入的复数组进行快速傅里叶变换（FFT）
输入参数：
*xin复数结构体组的首地址指针，struct型
*****************************************************************/
voidFFT（structcompx*xin）
{
intf,m,nv2,nm1,i,k,l,j=0;
structcompxu,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
}
{
intle,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）;//改变系数，进行下一个蝶形运算
}
}
}
}
/************************************************************
函数原型：
voidmain（）
函数功能：
测试FFT变换，演示函数使用方法
输入参数：
无
输出参数：
无
************************************************************/
voidmain（）
{
inti;
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;is[i].real=sqrt（s[i].real*s[i].real+s[i].imag*s[i].imag）;
while
（1）;
}
%////二、
%/////////////////////////////////////////////////////////////////////////////////////////////////////////////
%///////////////////////////////////////////////////////////////////////////////////////////////////////////
%////////////////////////////////MATLAB仿真信号源的源程序:
:
Clear;
Clc;
t=O:
O.01:
3;
yl=100*sin（pi/3*t）;
n=l;
fort=-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（'Nodenumber'）
T=input（'caiyangjiange'）
f=input（'frenquency'）
choise=input（'addzeroornot?
1/0'）
n=0:
T:
（N-1）*T;%采样点
k=0:
N-1;
x=sin（2*pi*f*n）;
ifchoise==1
e=input（'Inputthenumberofzeros!
'）
x=[xzeros（1,e）]
N=N+e;
else
end%加零
k=0:
N-1;%给k重新赋值，因为有可能出现加零状况
bianzhi=bi2de（fliplr（de2bi（k,length（de2bi（N-1）））））+1;%利用库函数进行变址运算
forl=1:
N
X（l）=x（bianzhi（l））;%将采样后的值按照变址运算后的顺序放入X矩阵中
end
d1=1;
form=1:
log2（N）
d2=d1;%做蝶形运算的两个数之间的距离
d1=d1*2;%同一级之下蝶形结之间的距离
W=1;%蝶形运算系数的初始值
dw=exp（-j*pi/d2）;%蝶形运算系数的增加量
fort=1:
d2%
fori=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
#defineN1000
typedefstruct
{
doublereal;
doubleimg;
}complex;
voidfft（）;/*快速傅里叶变换*/
voidifft（）;/*快速傅里叶逆变换*/
voidinitW（）;
voidchange（）;
voidadd（complex,complex,complex*）;/*复数加法*/
voidmul（complex,complex,complex*）;/*复数乘法*/
voidsub（complex,complex,complex*）;/*复数减法*/
voiddivi（complex,complex,complex*）;/*复数除法*/
voidoutput（）;/*输出结果*/
complexx[N],*W;/*输出序列的值*/
intsize_x=0;/*输入序列的长度，只限2的N次方*/
doublePI;
intmain（）
{
inti,method;
system（"cls"）;
PI=atan
（1）*4;
printf（"Pleaseinputthesizeofx:
\n"）;
/*输入序列的长度*/
scanf（"%d",&size_x）;
printf（"Pleaseinputthedatainx[N]:
（suchas:
56）\n"）;
/*输入序列对应的值*/
for（i=0;iscanf（"%lf%lf",&x[i].real,&x[i].img）;
initW（）;
/*选择FFT或逆FFT运算*/
printf（"UseFFT（0）orIFFT
（1）?
\n"）;
scanf（"%d",&method）;
if（method==0）
fft（）;
else
ifft（）;
output（）;
return0;
}
/*进行基-2FFT运算*/
voidfft（）
{
inti=0,j=0,k=0,l=0;
complexup,down,product;
change（）;
for（i=0;i（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;
}
}
}
}
voidifft（）
{
inti=0,j=0,k=0,l=size_x;
complexup,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（）;
}
voidinitW（）
{
inti;
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）;
}
}
voidchange（）
{
complextemp;
unsignedshorti=0,j=0,k=0;
doublet;
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;
}
}
}
voidoutput（）/*输出结果*/
{
inti;
printf（"Theresultareasfollows\n"）;
for（i=0;i{
printf（"%.4f",x[i].real）;
if（x[i].img>=0.0001）
printf（"+%.4fj\n",x[i].img）;
elseif（fabs（x[i].img）<0.0001）
printf（"\n"）;
else
printf（"%.4fj\n",x[i].img）;
}
}
voidadd（complexa,complexb,complex*c）
{
c->real=a.real+b.real;
c->img=a.img+b.img;
}
voidmul（complexa,complexb,complex*c）
{
c->real=a.real*b.real-a.img*b.img;
c->img=a.real*b.img+a.img*b.real;
}
voidsub（complexa,complexb,complex*c）
{
c->real=a.real-b.real;
c->img=a.img-b.img;
}
voiddivi（complexa,complexb,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*cos（2*pi*113.8.*t）+50*cos（2*pi*114.5.*t）;
[f, y] = zfft（x, 109, 115, fs）;
plot（f, y）;
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
图片效果：
升级会员 