用C实现FFTWord格式.docx
《用C实现FFTWord格式.docx》由会员分享，可在线阅读，更多相关《用C实现FFTWord格式.docx（16页珍藏版）》请在冰点文库上搜索。
xlabel（'
Time（s）'
ylabel（'
Amplitude（volt）'
title（'
Singnal'
holdon;
stem（t2,signal2,'
r'
axis（[0T*（N+ZeroNum）-11]）;
%作FFT变换，计算其幅值，归一化处理，并画出频谱。
Y=fft（signal2,N）;
Pyy=Y.*conj（Y）;
Pyy=（Pyy/sum（Pyy））*2;
f=0:
fs/（N-1）:
fs/2;
4
subplot（2,1,2）;
bar（f,Pyy（1:
N/2））;
Frequency（Hz）'
Amplitude'
Frequencycompnentsofsignal'
axis（[0fs/20ceil（max（Pyy））]）
gridon;
———————————————————————————————————————
这是一个傅里叶变化的子函数，你可以自己做主函数传递你这里的参数验证
//入口参数：
//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,返回逆傅立叶变换的虚部
il=1,l=0时，返回傅立叶变换的模
//il=1,l=1时，返回逆傅立叶变换的模
il=1,l=0时，返回傅立叶变换的辐角
//il=1,l=1时，返回逆傅立叶变换的辐角
voidkbfft（double*pr,double*pi,intn,intk,double*fr,double*fi,intl,intil）
{
intit,m,is,i,j,nv,l0;
doublep,q,s,vr,vi,poddr,poddi;
//排序
for（it=0;
it<
=n-1;
it++）
{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/（1.0*n）;
pr[1]=cos（p）;
pi[1]=-sin（p）;
if（l!
=0）pi[1]=-pi[1];
for（i=2;
{p=pr[i-1]*pr[1];
q=pi[i-1]*pi[1];
s=（pr[i-1]+pi[i-1]）*（pr[1]+pi[1]）;
pr[i]=p-q;
pi[i]=s-p-q;
=n-2;
it=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=m/2;
nv=2*nv;
=（m-1）*nv;
it=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=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]=fr[it+j]+poddr;
fi[it+j]=fi[it+j]+poddi;
=0）
{fr[i]=fr[i]/（1.0*n）;
fi[i]=fi[i]/（1.0*n）;
if（il!
{pr[i]=sqrt（fr[i]*fr[i]+fi[i]*fi[i]）;
pr[i]=（pr[i]/（n/2））;
//各次谐波幅值，其中pr[1]为基波幅值
if（fabs（fr[i]）<
0.000001*fabs（fi[i]））//fabs（）是取绝对值函数,浮点型的0在内存中并不是严格等于0,可以认为当一个浮点数离原点足够近时,也就是f>
0.00001&
&
f<
-0.00001,认为f是0
{if（（fi[i]*fr[i]）>
0）pi[i]=90.0;
elsepi[i]=-90.0;
else
pi[i]=atan（fi[i]/fr[i]）*360.0/6.283185306;
return;
}
FFTinCcode用C语言实现FFT
#include<
stdio.h>
math.h>
//此代码来源《数字信号处理C语言程序集》殷福亮、宋爱军，沈阳：
辽宁科学技术出版社，1997.7
//数组x存储时域序列的实部，数组y存储时域序列的虚部
//n代表N点FFT，sign=1为FFT，sign=-1为IFFT
voidFFT（doublex[],doubley[],intn,intsign）
{
inti,j,k,l,m,n1,n2;
doublec,c1,e,s,s1,t,tr,ti;
//Calculatei=log2N
for（j=1,i=1;
16;
m=i;
j=2*j;
if（j==n）
break;
//计算蝶形图的输入下标（码位倒读）
n1=n-1;
for（j=0,i=0;
n1;
if（i<
j）
{
tr=x[j];
ti=y[j];
x[j]=x[i];
y[j]=y[i];
x[i]=tr;
y[i]=ti;
k=n/2;
while（k<
（j+1））
j=j-k;
k=k/2;
j=j+k;
//计算每一级的输出，l为某一级，i为同一级的不同群，使用同一内存（即位运算）
n1=1;
for（l=1;
l<
=m;
l++）
n1=2*n1;
n2=n1/2;
e=3.1415926/n2;
c=1.0;
s=0.0;
c1=cos（e）;
s1=-sign*sin（e）;
for（j=0;
n2;
for（i=j;
n;
i+=n1）
k=i+n2;
tr=c*x[k]-s*y[k];
ti=c*y[k]+s*x[k];
x[k]=x[i]-tr;
y[k]=y[i]-ti;
x[i]=x[i]+tr;
y[i]=y[i]+ti;
t=c;
c=c*c1-s*s1;
s=t*s1+s*c1;
//如果是求IFFT，再除以N
if（sign==-1）
for（i=0;
x[i]/=n;
y[i]/=n;
intmain（）
doublex[4]={1,2,-1,3};
//郑君里《信号与系统》上的例子，下P145
doubley[4]={0};
inti;
4;
printf（"
%d%f%f\n"
i,x[i],y[i]）;
FFT（x,y,4,1）;
FFT（x,y,4,-1）;
return1;
《数字信号处理C语言程序集》殷福亮、宋爱军
pdf下载地址
以前信号与系统学得很浅，目前正在学习，代码只是看懂了个大概，对于是如何实现W的计算的没看懂，希望有谁看懂了解释一下，呵呵......
还有一点我没有看懂就是，变换出来的频域序列和频率是怎么对应的，好像是N/2=采样频率的一半，正好满足耐奎斯特，N/2以后和以前是对称的，不知道理解的正确与否，希望牛人指点....
C语言FFT和IFFT实现（2007-05-0718:
40:
15）转载标签：
基_2fft（频率抽取dif法）c语言算法程序分类：
学习资料
基-2FFT（频率抽取DIF法）算法程序
/*基2-DIF-FFT的C语言算法程序*/
/*输入:
序列点数、序列值*/
/*输出:
输入序列FFT变换后的数值及又把FFT结果反变换得到原输入序列（与输入序列相同），*/
/*还输出输入序列的IFFT结果（应与原序列相同）*/
#include"
conio.h"
math.h"
stdio.h"
#defineN64
#definePI3.1415926
#definew0（0.125*PI）
#defineCmul（a,b,c）a.x=b.x*c.x-b.y*c.y;
a.y=b.x*c.y+b.y*c.x;
/*利用结构体做复数乘法*/
#defineCequal（a,b）a.x=b.x;
a.y=b.y;
/*结构体a的值赋给b*/
#defineCadd（a,b,c）a.x=b.x+c.x;
a.y=b.y+c.y;
/*结构体实现复数的加减*/
#defineCsub（a,b,c）a.x=b.x-c.x;
a.y=b.y-c.y;
#defineWn（w,r）w.x=cos（2.0*PI*r/n）;
w.y=-sin（2.0*PI*r/n）;
/*旋转因子的分解*/
structcomp
floatx;
floaty;
};
voidmain（）
inti,j,nu2,nm1,n,m,le,le1,k,ip,z;
intflag,f,n1;
structcompa[N],t,t1,w,d;
floata_ipx,m1;
printf（"
\nThisprogramisaboutFFTbyDIFway."
\npleaseenterN:
"
scanf（"
%d"
&
n1）;
n=n1;
m1=log（n1）/log
（2）;
m=log（n1）/log
（2）;
/*级数*/
if（m!
=m1）n=pow（2,m+1）;
for（i=0;
i<
i++）{a[i].x=a[i].y=0.0;
}/*所有存储单元清零*/
\n"
i++）/*录入数据到存储单元，以结构体的形式分别存储实部和虚部*/
\npleaseenterdata（%d）_[Re]:
i）;
%f"
a[i].x）;
\npleaseenterdata（%d）_[Im]:
a[i].y）;
for（z=0;
z<
=2;
z++）/*z=0,时计算FFT；
z=1时计算新序列的IFFT；
z=2时计算原序列的IFFT*/
flag=-1;
/*flag控制r的增加*/
for（m=（log（n）/log
（2））;
m>
=1;
m--）/*此循环控制级数m*/
le=pow（2,m）;
/*le控制同一级中具有相同的碟形结的距离*/
flag++;
le1=le/2;
/*le1控制在同一级中有不同的个数*/
for（j=0;
j<
le1;
j++）/*此循环控制同一级中r的增加，实现一级不同的蝶形运算*/
for（i=j;
=（n-1）;
i+=le）/*此循环控制一级中具有相同的碟形运算*/
ip=i+le1;
Cequal（t,a[i]）;
Cequal（t1,a[ip]）;
f=（int）（i*pow（2,flag））%n;
Wn（w,f）;
Cadd（a[i],t,t1）;
/*计算*/
Csub（a[ip],t,t1）;
/*计算
a_ipx=a[ip].x;
*/
if（z==1）
w.y*=-1;
if（z==2）
a[ip].x=a[ip].x*w.x-a[ip].y*w.y;
a[ip].y=a_ipx*w.y+a[ip].y*w.x;
/*计算，反变换则取反*/
nu2=n/2;
/*用雷德算法将倒位序输出变为自然序输出*/
nm1=n-2;
j=0;
i=0;
while（i<
=nm1）
if（i<
j）
Cequal（d,a[j]）;
Cequal（a[j],a[i]）;
Cequal（a[i],d）;
k=nu2;
=j）
j=j-k;
k=k/2;
j=j+k;
i=i+1;
if（z==0）
\ninputnumfft:
\n\n"
elseif（z==1）
\nshengchengxinnumifft（input）:
）;
\ninputnumifft"
i++）
%7.4f"
a[i].x）;
if（a[i].y>
+%7.4fj\n"
a[i].y）;
-%7.4fj\n"
fabs（a[i].y））;
a[i].x=a[i].x/n;
/*IFFT需将结果除以n*/
a[i].y=a[i].y/n;
a[i].x/n）;
升级会员 