数字信号处理实验FFT快速傅里叶变换C语言.docx

1、数字信号处理实验FFT快速傅里叶变换C语言数字信号处理实验是按正常顺序排列在存储单元中，即按X(0,X(1,X(7的顺序排列，但是这时输入x(n却不是按自然顺序存储的，而是按x(0,x(4,x(7的顺序存入存储单元，所以我们要对输入的按正常顺序排列的数据进行变址存储，最后才能得到输出的有序的X(K。通过观察，可以发现，如果说输出数据是按原位序排列的话，那么输入数据是按倒位序排列的。即如果输入序列的序列号用二进制数,则到位序就为。我们需将输入的数据变为输出的倒位序存储，这里用雷德算法来实现。下面给出雷德算法。假如使用AI存的是顺序位序，而BJ存的是倒位序。例如 N = 8 的时候，倒位序顺序 二

2、进制表示 倒位序顺序0 0 000 0004 1 100 0012 2 010 010 6 3 110 0111 4 001 1005 5 101 1013 6 011 1107 7 111 111由上面的表可以看出，按自然顺序排列的二进制数，其下面一个数总是比其上面一个数大1，即下面一个数是上面一个数在最低位加1并向高位进位而得到的。而倒位序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位进位而得到。I、J都是从0开始，若已知某个倒位序J，要求下一个倒位序数，则应先判断J的最高位是否为0，这可与k=N/2相比较，因为N/2总是等于100.的。如果kJ，则J的最高位为0，只要把该位变

3、为1J与k=N/2相加即可），就得到下一个倒位序数；如果K=J，则J的最高位为1，可将最高位变为0J与k=N/2相减即可）。然后还需判断次高位，这可与k=N4相比较，若次高位为0，则需将它变为1加N4即可）其他位不变，既得到下一个倒位序数；若次高位是1，则需将它也变为0。然后再判断下一位。2. 复数运算因为每一个蝶形结构完成的迭代运算为算式中涉及到了复数的运算，而计算机是不能自己实现复数运算的，所以需要我们自己设计进行复数运算的程序。迭代运算式中，= cos2r/N）- jsin + jI(j, = R(K + jI(K ,而我们最后希望得到的DFT结果是复数的模，根据它的模来绘制频谱，所以这

4、里我们定义X(k=相关程序我们编译为：c.real=a.real*b.real-a.imag*b.imag。 c.imag=a.real*b.imag+a.imag*b.real。根据迭代运算的式子，我们可以将其分解为：R(K+jI(K=R(K+jI(K+ R(j + jI(j* cos2r/N）-jsin= R(K+ R(j cos sin (2r/N I(K=I(K-R(j sincos (2r/N 同理R(j= R(K- R(j cos sin (2r/N I(j=I(K+R(j sincos (2r/N 相关程序编译为：xinip.real=xini.real-t.real。 xini

5、p.imag=xini.imag-t.imag。 xini.real=xini.real+t.real。 xini.imag=xini.imag+t.imag。3. 节点距离运算当输入为倒位序，输出为正常顺序时，第m级运算，每个蝶形的两节点距离为。4. 旋转因子的计算这里主要解决的是迭代运算的每一项应乘的旋转因为中r应取多少的问题。第L级的2(L-1个碟形因子WPN 中的P，可表示为p = j*2(m-L，其中j = 0，1，2，.，个蝶形因子，第二层循环根据乘数进行控制，保证对于每一个旋转因子第三层循环要执行一次，这样，第三层循环在第二层循环控制下，每一级要进行2(L-1次循环计算。第三层：

6、因为第L级共有N/2L个蝶形结构，并且同一级内不同蝶型结构的旋转因子分布相同，当第二层循环确定某一旋转因子后，第三层循环要将本级中每个蝶型结构中具有这一旋转引自的蝶形计算一次，即第三层循环每执行完一次要进行N/2L个蝶形计算。所以，在每一级中，第三层循环完成N/2L个蝶形计算；第二层循环使得第三层循环进行 2(L-1次，因此，第二层循环完成时，共进行2(L-1 *N/2L=N/2个碟形计算。实质是：第二、第三层循环完成了第L级的计算。五、程序代码1.C语言：#include#include#include#define PI 3.14932384626433832795028841971 #d

7、efine FFT_N 128 /定义傅利叶变换的点数struct compx double real,imag。 /定义一个复数结构struct compx sFFT_N。 /从S1开始存放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 /FFT函数 int f,m,nv2,nm1,i

8、,k,l,j=0。 struct compx u,w,t。 nv2=FFT_N/2。 /变址运算，即把自然顺序变成倒位序，采用雷德算法 nm1=FFT_N-1。 for(i=0。i if(i /如果ij,即进行变址 t=xinj。 xinj=xini。 xini=t。 k=nv2。 /求j的下一个倒位序 while(k /如果k!=1。l+ /第一层循环，计算蝶形级数 。 for(m=1。m /第二层循环，控制蝶形级数 d=2。 /d蝶形结距离，即第m级蝶形的蝶形结构相距d点 d1=d/2。 /同一蝶形结中参加运算的两点的距离 u.real=1.0。 /u为蝶形结构运算系数，初始值为1 u.i

9、mag=0.0。 w.real=cos(PI/d1。 /w为系数商，即当前系数与前一个系数的商 w.imag=-sin(PI/d1。 for(j=0。j /第三层循环，进行旋转因子的计算完成蝶形计算。这里控制计算系数不同的蝶形结 for(i=j。i /控制计算系数相同蝶形结 ip=i+d1。 /i，ip分别表示参加蝶形运算的两个节点 t=EE(xinip,u。 /蝶形运算 xinip.real=xini.real-t.real。 xinip.imag=xini.imag-t.imag。 xini.real=xini.real+t.real。 xini.imag=xini.imag+t.imag

10、。 u=EE(u,w。 /改变系数，进行下一个蝶形运算 void main( int i。 for(i=0。i /给结构体赋值 si.real=sin(2*3.1493*i/FFT_N。 si.imag=0。 long start=clock(。 FFT(s。 for(i=0。i si.real=sqrt(si.real*si.real+si.imag*si.imag。long end=clock(。 printf(%.4fn,si.real。long t=end-start。printf(%dn,t。2.MATLAB：n=0:127。b=sin(2*pi*n/128。fft(b六、实验结果1.C语言2.MATLAB Columns 1 through 5 0.0000 -0.0000 -64.0000i -0.0000 - 0.0000i -0.0000