FFT算法实现实验报告
辛旸 PB10210006
实验目的
1、加深对快速傅里叶变换的理解。
2、掌握FFT 算法及其程序的编写。
3、掌握算法性能评测的方法。
实验内容
1.编写自己的FFT算法：
代码如下：
function [ X ] = Sampling( x,N )
%myFFT 实现FFT时域取样算法
% 输出：生成FFT序列X(k),输入：欲变换序列x（n），FFT变换长度N（可缺省）
（1） if ~exist('N','var'); %检查是否有变换长度N输入
（2） N=length(x); %若无，则令N等于序列长度
（3） end
（4） if N<length(x); %如果N小于序列长度，则对序列进行截短
（5） x=x(:,1:N);
（6） else
（7） x=[x,zeros(1,N-length(x))]; %如果N大于序列长度，对序列补零进行延长
（8） end;
（9） for i=1:1:length(x)/2+1; %判断N是2的多少次方
（10） if 2^i>=length(x); %若N不是2的整数幂
（11） N=2^i; %增大N为2的整数幂
（12） break;
（13） end
（14） end
（15） x=[x,zeros(1,N-length(x))]; %确保要变换的序列长度为2^i
（16） k1=zeros(1,N);
（17） X=zeros(1,N);
（18） w=zeros(1,N);
（19） for m=0:1:N-1; %确定反序序列k1和正序序列k的关系
（20） k=m;
（21） for n=i-1:-1:0; %从高位开始依次将各位移至反序位
（22） k1(m+1)=k1(m+1)+fix(k/(2^n))*(2^(i-1-n));
（23） k=rem(k,2^n);
（24） end;
（25） end
（26） for l=1:1:N;
（27） X(k1(l)+1)=x(l); %生成反序序列
（28） w(l)=exp(-1i*2*pi/N*(l-1)); %生成旋转因子
（29） end
（30） for l=0:1:i-1; %控制FFT运算级数
（31） for m=1:1:N; %每一级中有N/2个蝶形运算
（32） if rem((m-1),2^(l+1))<2^l; %找到蝶形运算的上半部分
（33） b=X(m)+X(m+2^l)*w(2^(i-1-l)*rem((m-1),2^l)+1); %将结果暂存至b
（34） X(m+2^l)=a(m)-X(m+2^l)*w(2^(i-1-l)*rem((m-1),2^l)+1);
（35） X(m)=b; %实现原位运算
（36） end
（37） end
（38） end
2.选择实验1中的典型信号序列验证算法的有效性：
为方便比较两个算法，编写了myCompare函数计算两种算法的运行时间，并绘制频谱曲线
代码如下：
function [ t1,t2,e ] = myCompare( x,N )
%myCompare函数：比较自己编写的算法与系统自带算法的差异
% 输入：与变换信号序列x（n）和欲变换长度N
% 输出：自己编写的函数的执行时间t1，系统自带函数的执行时间t2，两者计算序列的差异平方和e tic;
X1=myFFT(x,N);
t1=toc;
tic;
X2=fft(x,N);
t2=toc;
subplot(1,2,1);plot(abs(X1));xlabel('k');ylabel('X(k');title('ÓÃ×Ô¼º±àдµÄº¯ÊýµÃµ½µÄ±ä»»ÐòÁÐƵÆ×');
subplot(1,2,2);plot(abs(X2));xlabel('k');ylabel('X(k');title('ϵͳ×Ô´øFFTº¯ÊýµÃµ½µÄ±ä»»ÐòÁÐƵÆ×');
e=sum((X1-X2).^2);
end
对理想采样信号A=444.128，α=50*2^(1/2)*π，Ω=50*2^(1/2)*π，T=1/1000，序列长度50，用自己编写的FFT算法和系统自带算法做64点FFT变换后绘制频域序列，如下：
对高斯序列，p=8，q=8，序列长度16，用自己编写的FFT算法和系统自带算法做16点FFT变换后绘制频域序列，如下：
对衰减正弦序列α=0.01，f=0.05，序列长度100，用自己编写的FFT算法和系统自带算法做128点FFT变换后绘制频域序列，如下：
由以上结果可知，自己编写的算法运行结果与系统自带算法一致，且可以对信号进行截断或补零后再做变换。
3.对所编制FFT算法进行性能评估：
与自己编写的DFT算法进行性能比较：
对N点序列进行DFT变换需要N²/2次复乘，而对N点序列做基2-FFT只需N/2*log2（N）次复乘，因此运算量减少了很多，且随着序列长度增加，运算量差异变大。
与系统自带FFT算法进行性能比较：
由于系统自带FFT函数用C语言实现，无法查看源代码，只知道效率更高，而且在计算任意点的DFT（不指定变换长度N）时，系统自带函数无需采取补零操作，而自己编写的函数会先补零再变换，改变了频域取样密度，会得到与系统自带函数不同的结果。
比较自己编写的DFT算法、自己编写的FFT算法和系统自带算法三者运算时间，得到下表：
由此绘制曲线如图：
由比较结果可知，虽然运算时间曲线和理想曲线不完全相同，但三种算法的相对运算时间与理论预期一致。
实验总结及个人结论：
从对实验的比较结果中可知，自己编写的FFT算法有效且比自己编写的DFT算法快很多，但却始终比系统自带的算法慢，一方面是因为系统自带算法实现效率高，另一方面，在观察了自己编写的算法中各步执行时间后，发现函数用在生成反序序列的时间和实际进行FFT运算的时间相当，相当于多用了一倍的时间。另外，对任意长序列，自己编写的FFT算法只能补零后计算，而不能像系统自带函数一样算出实际的N点DFT。
第二篇：基于DSP的FFT算法实现
基于DSP的FFT算法实现
1、 FFT的原理
快速傅氏变换（FFT）是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。
设x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2（N/2）2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。
数字信号处理器(DSP)是一种可编程的高性能处理器，近年来发展很快．它不仅适用于数字信号处理，而且在图像处理、语音处理、通信等领域得到了广泛的应用．通用的微处理器在运算速度上很难适应信号实时处理的要求．联沪处理器中集成有高速的乘法器硬件，能快速地进行大量数据的乘法和加法运算。快速傅里叶变换(FFT)的出现使得DFr在实际应用中得到了广泛的应用．
2、 基于DSP的FFT算法实现
用C语言实现FFT算法
/*****************fft programe*********************/
#include "typedef.h"
#include "math.h"
struct compx EE(struct compx b1,struct compx b2)
{
struct compx b3 ;
b3.real=b1.real*b2.real-b1.imag*b2.imag ;
b3.imag=b1.real*b2.imag+b1.imag*b2.real ;
return(b3);
}
void FFT(struct compx*xin,int N)
{
int f,m,nv2,nm1,i,k,j=1,l ;
/*int f,m,nv2,nm1,i,k,j=N/2,l;*/
struct compx v,w,t ;
nv2=N/2 ;
f=N ;
for(m=1;(f=f/2)!=1;m++)
{
;
}
nm1=N-1 ;
/*变址运算*/
for(i=1;i<=nm1;i++)
{
if(i<j)
{
t=xin[j];
xin[j]=xin[i];
xin[i]=t ;
}
k=nv2 ;
while(k<j)
{
j=j-k ;
k=k/2 ;
}
j=j+k ;
}
{
int le,lei,ip ;
float pi ;
for(l=1;l<=m;l++)
{
le=pow(2,l);
// 这里用的是L而不是1
lei=le/2 ;
pi=3.14159 ;
v.real=1.0 ;
v.imag=0.0 ;
w.real=cos(pi/lei);
w.imag=-sin(pi/lei);
for(j=1;j<=lei;j++)
{
/*double p=pow(2,m-l)*j;
double ps=2*pi/N*p;
w.real=cos(ps);
w.imag=-sin(ps);*/
for(i=j;i<=N;i=i+le)
{
/* w.real=cos(ps);
w.imag=-sin(ps);*/
ip=i+lei ;
t=EE(xin[ip],v);
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 ;
}
v=EE(v,w);
}
}
}
return ;
}
/*****************main programe********************/
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include "typedef.h"
float result[257];
struct compx s[257];
int Num=256 ;
const float pp=3.14159 ;
main()
{
int i=1 ;
for(;i<0x101;i++)
{
s[i].real=sin(pp*i/32);
s[i].imag=0 ;
}
FFT(s,Num);
for(i=1;i<0x101;i++)
{
result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));
}
}
3、ICETEK-F2812-A的实验板调试步骤
1．实验准备
(1)连接实验设备：
(2)准备信号源进行AD 输入。
①取出2 根实验箱附带的信号线(如右图，两端均为单声道语音插头)。
②用1 根信号线连接实验箱底板上信号源I模块(图10-1 中单实线框
中部分)的“波形输出”插座(图10-1中的3 或4)和“A/D 输入”模块(图10-1中虚线框
中部分)的“ADCIN0”插座(图10-1 中的A)，注意插头要插牢、到底。这样，信号源I
的输出波形即可送到ICETEK-F2812-A评估板的AD 输入通道0。
③用1 根信号线连接实验箱底板上信号源II模块(图10-1中双实线框中部分)的“波形输
出”插座(图10-1 中的c或d)和“A/D 输入”模块的“ADCIN1”插座(图10-1中的B)，
注意插头要插牢、到底。这样，信号源II的输出波形即可送到ICETEK-F2812-A评估
板的AD 输入通道1。
④设置信号源I：
-调整拨动开关“频率选择”(图10-1 中的5)拨到“100Hz－1KHz”档(图10-1中10)。
-将“频率微调”(图10-1 中的6)顺时针调到头(最大)。
-调整拨动开关“波形选择”(图10-1 中的7)拨到“三角波”档(图10-1 中的11)。
-将“幅值微调”（图10-1 中的8）顺时针调到头（最大）。
⑤设置信号源II：
-调整拨动开关“频率选择”(图10-1 中的e)拨到“100Hz－1KHz”档（图10-1 中的j）。
-将“频率微调”（图10-1 中的f）顺时针调到头（最大）。
-调整拨动开关“波形选择”（图10-1 中的g）拨到“正弦波”档（图10-1 中的k）。
-将“幅值微调”（图10-1 中的h）顺时针调到头（最大）。
⑥将两个信号源的电源开关(图10-1 中的2和b)拨到“开”的位置。
2．设置Code Composer Studio 2.21在硬件仿真（Emulator）方式下运行
请参看本书第一部分、四、2。
3．启动Code Composer Studio 2.21
请参看本书第一部分、五、2。
选择菜单Debug→Reset CPU。
4．打开工程文件
-工程目录：C:\ICETEK-F2812-A-EDUlab\DSP281x_examples\Lab0305-AD\ ADC.pjt。
-在项目浏览器中，双击adc.c，打开adc.c 文件，浏览该文件的内容，理解各语句作用。
5．编译、下载程序。
6．打开观察窗口
-选择菜单“View”、“Graph”、“Time/Frequency…”做如下设置，然后单击“OK”按钮；
-选择菜单“View”、“Graph”、“Time/Frequency…”做如下设置(图10-3)，然后单击“OK”
按钮；
-在弹出的图形窗口中单击鼠标右键，选择“Clear Display”。
通过设置，我们打开了两个图形窗口观察两个通道模数转换的结果。
7． 设置信号源
由于模数输入信号未经任何转换就进入DSP，所以必须保证输入的模拟信号的幅度在
0-3V之间。必须用示波器检测信号范围，保证最小值0V最大值3 V，否则容易损坏DSP
芯片的模数采集模块。
8．运行程序观察结果
-单击“Debug”菜单，“Run”项，运行程序；
-停止运行，观察“ADCIN0”、“ADCIN1”窗口中的图形显示；
-适当改变信号源，按F5 健再次运行，停止后观察图形窗口中的显示。
注意：输入信号的频率不能大于10KHz，否则会引起混叠失真，而无法观察到波形，如
果有兴趣，可以试着做一下，观察采样失真后的图形。
9．选择菜单File→workspace→save workspacs As…，输入文件名SY.wks 。
10．退出CCS
4、DSP板调试结果之波形及波形分析
DSP板有多种调试方法，下面使用的一种为观测其示波器波形的方法，通过波形，我们可以看出此板子的性能。使用通用定时器Timer1/2/3/4产生PWM，选择连续计数模式可以产生如下图所示的非对称PWM波形
选择连续增/减计数模式可以产生中心或对称PWM波形，如下图所示
（2）同样，采用连续增计数模式可以产生一对带有死区的互补的非对称PWM波形
采用连续增/减计数模式可以产生一对带有死区的互补的对称PWM波形
实现的方法如下：
（1） 采用通用定时器Timer1和Timer2产生两路PWM波形；
（2） 为了产生对称波形，使两个定时器都工作于连续增/减计数模式；
（3） 从上图可以看出，S1的上升延和S2的上升延始终相差半个Ts，及半个周期，为了实现相移，可以让T1先开始计时工作，当T1到达第一个周期中断的时候打开T2，让T2也开始工作，同时需要去使能T1中断，或者通过置标志位等方法使得以后T1周期中断的时候不会再去打开T2定时器。这样就可以使的T1和T2输出的波形满足上图的要求，即既是互补，又是导通时间对称的PWM波形，只要占空比不足50％，就相当于留有一段死区时间。可以参看下面的示意图。