C语言Matlab实现FFT几种编程实例Word格式.docx
《C语言Matlab实现FFT几种编程实例Word格式.docx》由会员分享，可在线阅读，更多相关《C语言Matlab实现FFT几种编程实例Word格式.docx（14页珍藏版）》请在冰豆网上搜索。
nm1=FFT_N-1;
for（i=0;
i<
nm1;
i++）
{
if（i<
j）//如果i<
j,即进行变址
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;
=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变换，演示函数使用方法
无
************************************************************/
{
inti;
FFT_N;
i++）//给结构体赋值
s[i].real=sin（2*3.141592653589793*i/FFT_N）;
//实部为正弦波FFT_N点采样，赋值为1
s[i].imag=0;
//虚部为0
FFT（s）;
//进行快速福利叶变换
i++）//求变换后结果的模值，存入复数的实部部分
s[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:
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重新赋值，因为有可能出现加零状况
bianzhi=bi2de（fliplr（de2bi（k,length（de2bi（N-1）））））+1;
%利用库函数进行变址运算
forl=1:
N
X（l）=x（bianzhi（l））;
%将采样后的值按照变址运算后的顺序放入X矩阵中
d1=1;
form=1:
log2（N）
d2=d1;
%做蝶形运算的两个数之间的距离
d1=d1*2;
%同一级之下蝶形结之间的距离
W=1;
%蝶形运算系数的初始值
dw=exp（-j*pi/d2）;
%蝶形运算系数的增加量
fort=1:
d2%
fori=t:
d1:
i1=i+d2;
if（i1>
N）break;
%判断是否超出范围
else
p=X（i1）*W;
X（i1）=X（i）-p;
X（i）=X（i）+p;
%蝶形运算
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<
stdio.h>
stdlib.h>
#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）;
Pleaseinputthedatainx[N]:
（suchas:
56）\n"
/*输入序列对应的值*/
size_x;
i++）
%lf%lf"
x[i].real,&
x[i].img）;
initW（）;
/*选择FFT或逆FFT运算*/
UseFFT（0）orIFFT
（1）?
method）;
if（method==0）
fft（）;
else
ifft（）;
output（）;
return0;
}
/*进行基-2FFT运算*/
voidfft（）
inti=0,j=0,k=0,l=0;
complexup,down,product;
change（）;
log（size_x）/log
（2）;
l=1<
i;
j+=2*l）
for（k=0;
k<
l;
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;
（int）（log（size_x）/log
（2））;
i++）/*蝶形运算*/
l/=2;
j+=2*l）
k++）
add（x[j+k],x[j+k+l],&
up.real/=2;
up.img/=2;
sub（x[j+k],x[j+k+l],&
down.real/=2;
down.img/=2;
divi（down,W[size_x*k/2/l],&
voidinitW（）
W=（complex*）malloc（sizeof（complex）*size_x）;
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;
k=i;
j=0;
t=（log（size_x）/log
（2））;
while（（t--）>
0）
j=j<
1;
j|=（k&
1）;
k=k>
>
if（j>
i）
temp=x[i];
x[i]=x[j];
x[j]=temp;
voidoutput（）/*输出结果*/
Theresultareasfollows\n"
%.4f"
x[i].real）;
if（x[i].img>
=0.0001）
+%.4fj\n"
x[i].img）;
elseif（fabs（x[i].img）<
0.0001）
%.4fj\n"
voidadd（complexa,complexb,complex*c）
c->
real=a.real+b.real;
img=a.img+b.img;
voidmul（complexa,complexb,complex*c）
real=a.real*b.real-a.img*b.img;
img=a.real*b.img+a.img*b.real;
voidsub（complexa,complexb,complex*c）
real=a.real-b.real;
img=a.img-b.img;
voiddivi（complexa,complexb,complex*c）
real=（a.real*b.real+a.img*b.img）/（
b.real*b.real+b.img*b.img）;
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:
b
2*pi*f0.*r
./
fs;
x1
x
.*
exp（-1j
b）;
%移频
bw
fa
-
fi;
B
fir1（32,
fs）;
%滤波
截止频率为0.5bw
x2
filter（B,
1,
x1）;
c
x2（1:
floor（fs/bw）:
N）;
N1
length（c）;
f
linspace（fi,
N1）;
y
abs（fft（c））
*
circshift（y,
[0,
floor（N1/2）]）;
%将负半轴的幅值移过来
end
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
%上边四程序测试应用实例：
fs
2048;
T
100;
t
1/fs:
T;
30
cos（2*pi*110.*t）
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）;
109,
115,
plot（f,
y）;
图片效果：
升级会员 