数字信号处理实验FFT快速傅里叶变换C语言Word格式文档下载.docx

数字信号处理实验FFT快速傅里叶变换C语言Word格式文档下载.docx

《数字信号处理实验FFT快速傅里叶变换C语言Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验FFT快速傅里叶变换C语言Word格式文档下载.docx（12页珍藏版）》请在冰豆网上搜索。

却不是按自然顺序存储的，而是按x（0>

x（4>

…,x（7>

的顺序存入存储单元，所以我们要对输入的按正常顺序排列的数据进行变址存储，最后才能得到输出的有序的X（K>

。

通过观察，可以发现，如果说输出数据是按原位序排列的话，那么输入数据是按倒位序排列的。

即如果输入序列的序列号用二进制数,则到位序就为。

我们需将输入的数据变为输出的倒位序存储，这里用雷德算法来实现。

下面给出雷德算法。

假如使用A[I]存的是顺序位序，而B[J]存的是倒位序。

例如 

N=8的时候，倒位序顺序 

二进制表示 

倒位序顺序

00 

000 

000

41 

100 

001

22 

010 

63 

110 

011

14 

001 

100

55 

101 

101

36 

011 

110

77 

111 

111

由上面的表可以看出，按自然顺序排列的二进制数，其下面一个数总是比其上面一个数大1，即下面一个数是上面一个数在最低位加1并向高位进位而得到的。

而倒位序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位进位而得到。

I、J都是从0开始，若已知某个倒位序J，要求下一个倒位序数，则应先判断J的最高位是否为0，这可与k=N/2相比较，因为N/2总是等于100..的。

如果k>

J，则J的最高位为0，只要把该位变为1<

J与k=N/2相加即可），就得到下一个倒位序数；

如果K<

=J，则J的最高位为1，可将最高位变为0<

J与k=N/2相减即可）。

然后还需判断次高位，这可与k=N\4相比较，若次高位为0，则需将它变为1<

加N\4即可）其他位不变，既得到下一个倒位序数；

若次高位是1，则需将它也变为0。

然后再判断下一位

2.复数运算

因为每一个蝶形结构完成的迭代运算为

算式中涉及到了复数的运算，而计算机是不能自己实现复数运算的，所以需要我们自己设计进行复数运算的程序。

迭代运算式中，=cos<

2πr/N）-jsin<

2πr/N）我们设=R（j>

+jI（j>

=R（K>

+jI（K>

而我们最后希望得到的DFT结果是复数的模，根据它的模来绘制频谱，所以这里我们定义

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>

+[R（j>

]*[cos<

2πr/N）-jsin<

2πr/N）]。

继续分解得到下列两式:

R（K>

=R（K>

+R（j>

cos<

2πr/N）+I（j>

sin（2πr/N>

I（K>

=I（K>

-R（j>

sin<

2πr/N）+I（j>

cos（2πr/N>

同理

R（j>

-R（j>

2πr/N）-I（j>

I（j>

+R（j>

2πr/N）I（j>

相关程序编译为：

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。

3.节点距离运算

当输入为倒位序，输出为正常顺序时，第m级运算，每个蝶形的两节点距离为。

4.旋转因子的计算

这里主要解决的是迭代运算的每一项应乘的旋转因为中r应取多少的问题。

第L级的2^（L-1>

个碟形因子WPN中的P，可表示为p=j*2^（m-L>

，其中j=0，1，2，...，<

2L-1-1）。

将上述方法整理一下，我们可以得到一个完整的程序的思路。

主要分为三个模块，复数计算DFT模函数，FFT函数，主函数。

其中复数计算DFT模函数的算法我们已经在前面介绍了。

主函数是编译整个程序进行的顺序，即我们先设定一个时域函数，从中抽取N点的值，然后进行FFT运算，然后求模长，最后将对应的频域的值输出。

这里比较复杂的是FFT函数。

下面进行解释。

这里我们共设三个循环来实现FFT这一函数的功能。

第一层：

首先我们知道，一个N=2^m点DFT，要进行m级计算,所以第一层循环是对运算的级数进行控制。

　　第二层：

因为第L级有2^（L-1>

个蝶形因子，第二层循环根据乘数进行控制，保证对于每一个旋转因子第三层循环要执行一次，这样，第三层循环在第二层循环控制下，每一级要进行2^（L-1>

次循环计算。

　　第三层：

因为第L级共有N/2^L个蝶形结构，并且同一级内不同蝶型结构的旋转因子分布相同，当第二层循环确定某一旋转因子后，第三层循环要将本级中每个蝶型结构中具有这一旋转引自的蝶形计算一次，即第三层循环每执行完一次要进行N/2^L个蝶形计算。

　　所以，在每一级中，第三层循环完成N/2^L个蝶形计算；

第二层循环使得第三层循环进行2^（L-1>

次，因此，第二层循环完成时，共进行2^（L-1>

*N/2^L=N/2个碟形计算。

实质是：

第二、第三层循环完成了第L级的计算。

五、程序代码

1.C语言：

#include<

stdio.h>

math.h>

time.h>

#definePI3.14932384626433832795028841971

#defineFFT_N128//定义傅利叶变换的点数

structcompx{doublereal,imag。

}。

//定义一个复数结构

structcompxs[FFT_N]。

//从S[1]开始存放

structcompxEE（structcompxa,structcompxb>

//求复数的模长函数

{

structcompxc。

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

return（c>

}

voidFFT（structcompx*xin>

//FFT函数

intf,m,nv2,nm1,i,k,l,j=0。

structcompxu,w,t。

nv2=FFT_N/2。

//变址运算，即把自然顺序变成倒位序，采用雷德算法

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。

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

j=j+k。

//把0改为1

intd,d1,ip。

f=FFT_N。

for（l=1。

（f=f/2>

!

=1。

l++>

//第一层循环，计算蝶形级数

。

for（m=1。

m<

=l。

m++>

//第二层循环，控制蝶形级数

{

d=2<

（m-1>

//d蝶形结距离，即第m级蝶形的蝶形结构相距d点

d1=d/2。

//同一蝶形结中参加运算的两点的距离

u.real=1.0。

//u为蝶形结构运算系数，初始值为1

u.imag=0.0。

w.real=cos（PI/d1>

//w为系数商，即当前系数与前一个系数的商

w.imag=-sin（PI/d1>

for（j=0。

j<

=d1-1。

j++>

//第三层循环，进行旋转因子的计算完成蝶形计算。

这里控制计算系数不同的蝶形结

for（i=j。

=FFT_N-1。

i=i+d>

//控制计算系数相同蝶形结

ip=i+d1。

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

t=EE（xin[ip],u>

//蝶形运算

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

u=EE（u,w>

//改变系数，进行下一个蝶形运算

voidmain（>

{

inti。

FFT_N。

//给结构体赋值

s[i].real=sin（2*3.1493*i/FFT_N>

s[i].imag=0。

longstart=clock（>

FFT（s>

{s[i].real=sqrt（s[i].real*s[i].real+s[i].imag*s[i].imag>

longend=clock（>

printf（"

%.4f\n"

s[i].real>

longt=end-start。

printf（"

%d\n"

t>

2.MATLAB：

n=0:

127。

b=sin（2*pi*n/128>

fft（b>

六、实验结果

1.C语言

2.MATLAB

Columns1through5

0.0000-0.0000-64.0000i-0.0000-0.0000i-0.0000