FFT算法
/****************************************************************************/ /* */ /* 快速傅里叶变换/ 快速傅里叶逆变换测试*/
/* */ /* 2014年04月20日*/ /* */ /****************************************************************************/ #include
#include
#include "mathlib.h" // DSP 数学函数库
#include "dsplib.h" // DSP 函数库
/****************************************************************************/ /* */ /* 宏定义*/ /* */ /****************************************************************************/ // 软件断点
#define SW_BREAKPOINT asm(" SWBP 0 ");
// 快速傅里叶变换
// π及浮点数极小值
#define PI 3.14159
#define F_TOL (1e-06)
/****************************************************************************/ /* */ /* 全局变量*/ /* */ /****************************************************************************/ // 快速傅里叶变换测试
// 测试快速傅里叶变换点数
// 注意:TI DSP库最大支持一次性计算128K 个点的FFT
#define Tn 1024
// 采样频率
#define Fs 1000.0
// 信号
float Input[2*Tn+4];
// FFT 输入信号
#pragma DATA_ALIGN(CFFT_In, 8);
float CFFT_In[2*Tn+4];
// FFT 输入信号副本
float CFFT_InOrig[2*Tn+4];
// FFT 输出
#pragma DATA_ALIGN(CFFT_Out, 8);
float CFFT_Out[2*Tn+4];
// IFFT 输出
#pragma DATA_ALIGN(CFFT_InvOut, 8);
float CFFT_InvOut[2*Tn+4];
// 中间运算临时变量
float CTemp[2*Tn+4];
// 存储旋转因子
float Cw[2*Tn];
// 模
float Cmo[Tn+2];
// 二进制位翻转
#pragma DATA_ALIGN (brev, 8);
unsigned char brev[64]=
{
0x0, 0x20, 0x10, 0x30, 0x8, 0x28, 0x18, 0x38,
0x4, 0x24, 0x14, 0x34, 0xc, 0x2c, 0x1c, 0x3c,
0x2, 0x22, 0x12, 0x32, 0xa, 0x2a, 0x1a, 0x3a,
0x6, 0x26, 0x16, 0x36, 0xe, 0x2e, 0x1e, 0x3e,
0x1, 0x21, 0x11, 0x31, 0x9, 0x29, 0x19, 0x39,
0x5, 0x25, 0x15, 0x35, 0xd, 0x2d, 0x1d, 0x3d,
0x3, 0x23, 0x13, 0x33, 0xb, 0x2b, 0x1b, 0x3b,
0x7, 0x27, 0x17, 0x37, 0xf, 0x2f, 0x1f, 0x3f
};
/****************************************************************************/ /* */ /* 函数声明*/ /* */ /****************************************************************************/ // 产生旋转因子
void tw_gen(float *w, int n);
// FFT 测试
void FFTTest();
/****************************************************************************/ /* */ /* 主函数*/ /* */ /****************************************************************************/ int main(void)
{
// FFT 测试
FFTTest();
// 断点
SW_BREAKPOINT;
}
/****************************************************************************/ /* */ /* 快速傅里叶变换测试*/ /* */ /****************************************************************************/ // 产生旋转因子
void tw_gen(float *w, int n)
{
int i,j,k;
double x_t,y_t,theta1,theta2,theta3;
for(j=1,k=0;j<=n>>2;j=j<<2)
{
for(i=0;i
{
theta1=2*PI*i/n;
x_t=cos(theta1);
y_t=sin(theta1);
w[k]=(float)x_t;
w[k+1]=(float)y_t;
theta2=4*PI*i/n;
x_t=cos(theta2);
y_t=sin(theta2);
w[k+2]=(float)x_t;
w[k+3]=(float)y_t;
theta3=6*PI*i/n;
x_t=cos(theta3);
y_t=sin(theta3);
w[k+4]=(float)x_t;
w[k+5]=(float)y_t;
k+=6;
}
}
}
// 快速傅里叶变换
void FFTTest(void)
{
// 产生待测试信号
unsigned int i;
for (i=0;i Input[i]=5*sin(2*PI*150*(i/Fs))+15*sin(2*PI*350*(i/Fs)); // 确定快速傅里叶变换基 unsigned char rad; if(Tn==16 || Tn==64 || Tn==256 || Tn==1024 || Tn==4096 || Tn==16384 || Tn==65536) rad=4; else if(Tn==8 || Tn==32 || Tn==128 || Tn==512 || Tn==2048 || Tn==8192 || Tn==32768) rad=2; else { printf ("不支持计算%d 点快速傅里叶变换！\n",Tn); return; } // 复数FFT for (i=0;i<2*Tn;i++) CFFT_In[i]=0.0; for (i=0;i { CFFT_In[2*i]=Input[i]; // 实部 CFFT_In[2*i+1]=0; // 虚部为0 } // 保留一份输入信号副本 memcpy(CFFT_InOrig,CFFT_In,2*Tn*sizeof(float)); // 产生旋转因子 tw_gen(Cw,Tn); // FFT 计算 DSPF_sp_fftSPxSP(Tn,CFFT_In,Cw,CFFT_Out,brev,rad,0,Tn); // 计算振幅 for(i=0;i Cmo[i]=0.0; for(i=0;i { Cmo[i]=sqrtsp(CFFT_Out[2*i]*CFFT_Out[2*i]+CFFT_Out[2*i+1]*CFFT_Out[2*i+1]); Cmo[i]=Cmo[i]*2/Tn; } // 保留一份FFT 结果副本 memcpy(CTemp,CFFT_Out,2*Tn*sizeof(float)); // IFFT 计算 DSPF_sp_ifftSPxSP(Tn,CFFT_Out,Cw,CFFT_InvOut,brev,rad,0,Tn); // 恢复FFT 结果 memcpy(CFFT_Out,CTemp,2*Tn*sizeof(float)); printf("\n复数FFT 测试结果:"); unsigned char Flag; for(i=0;i if(abs(CFFT_InOrig[i]-CFFT_InvOut[i])>F_TOL) Flag=1; if(Flag==1) printf ("失败！\n"); else printf ("成功！\n"); } 2.3 快速傅立叶变换问题 1) 问题背景 在数值电路的传输中，为了避免信号干扰，需要把一个连续信号 x(t)先通过取样离散化为一列数值脉冲信号x(0), x(1), …… ,然后再通过编码送到传输电路中。如果取样间隔很小，而连续信号的时间段又很长，则所得到的数值脉冲序列将非常庞大。因此，传输这个编码信号就需要长时间的占用传输电路，相应地也需要付出昂贵的电路费用。 那么能否经过适当处理是使上述的数值脉冲序列变短，而同时又不会丧失有用的信息？的经过研究，人们发现，如果对上述数值脉冲序列作如下的变换处理： ∑-=--=-== 1 /21 ,1,...,1,0,)()(N k N nki i N n e k x n X π (1) 则所得到的新序列X(0), X(1) , ……将非常有序，其值比较大的点往往集中在某一很狭窄的序列段内，这将非常有利于编码和存储，从而达到压缩信息的目的。 公式（1）就是所谓的离散傅立叶变换，简称DFT 。现在我们来分析一下计算DFT 所需要的工作量。如果我们不考虑公式（7.1）中指数项的运算，那么计算其每一个点X (n) 需要N 次复数乘法和N-1次的复数加法。显然当N 很大时，这个工作量也非常巨大。正是由于这个原因，使得DFT 的应用范围在过去很长的时间里受到了严格的限制。注意到公式（1）是非常有规律性的，那么能否利用这种规律性来降低DFT 的计算时间？ 1965年，凯莱和塔柯的提出了一种用于计算DFT 的数学方法，大大减少了DFT 的计算时间，同时又特别适用于硬件处理，这就是所谓的快速傅里叶变换，简称FFT 。鉴于DFT 的数据结构可以通过傅立叶变换的离散化获得，亦可通过三角插值得到，而本质上又同连续傅里叶分析有着极为密切的关系。下面我们从傅立叶级数级数和傅立叶积分入手，导出DFT 结构的来源和FFT 的工作原理。 2） 傅立叶变换 如果x(t)是定义在整个实轴上的实值或复值函数，则其傅立叶变换可由下式给出： ? ∞ ∞ ---== 1 ,)()(/2i dt e t x f X T nift （2） FFT是离散傅立叶变换的快速算法，可以将一个信号变换 到频域。有些信号在时域上是很难看出什么特征的，但是如 果变换到频域之后，就很容易看出特征了。这就是很多信号 分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱 提取出来，这在频谱分析方面也是经常用的。 虽然很多人都知道FFT是什么，可以用来做什么，怎么去 做，但是却不知道FFT之后的结果是什意思、如何决定要使用 多少点来做FFT。 现在圈圈就根据实际经验来说说FFT结果的具体物理意义。 一个模拟信号，经过ADC采样之后，就变成了数字信号。采样 定理告诉我们，采样频率要大于信号频率的两倍，这些我就 不在此罗嗦了。 采样得到的数字信号，就可以做FFT变换了。N个采样点， 经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT 运算，通常N取2的整数次方。 假设采样频率为Fs，信号频率F，采样点数为N。那么FFT 之后结果就是一个为N点的复数。每一个点就对应着一个频率 点。这个点的模值，就是该频率值下的幅度特性。具体跟原始 信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT 的结果的每个点（除了第一个点直流分量之外）的模值就是A 的N/2倍。而第一个点就是直流分量，它的模值就是直流分量 的N倍。而每个点的相位呢，就是在该频率下的信号的相位。 第一个点表示直流分量（即0Hz），而最后一个点N的再下一个 点（实际上这个点是不存在的，这里是假设的第N+1个点，也 可以看做是将第一个点分做两半分，另一半移到最后）则表示 采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率 依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。 由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果 采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时 间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率 第四章 快速傅里叶变换 有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换(FFT). 1965年，Cooley 和Tukey 提出了计算离散傅里叶变换（DFT ）的快速算法，将DFT 的运算量减少了几个数量级。从此，对快速傅里叶变换（FFT ）算法的研究便不断深入，数字信号处理这门新兴学科也随FFT 的出现和发展而迅速发展。根据对序列分解与选取方法的不同而产生了FFT 的多种算法，基本算法是基２DIT 和基２DIF 。FFT 在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。 快速傅里叶变换（FFT ）是计算离散傅里叶变换（DFT ）的快速算法。 DFT 的定义式为 )(k X =)()(1 0k R W n x N N n kn N ∑-= 在所有复指数值kn N W 的值全部已算好的情况下，要计算一个)(k X 需要N 次复数乘法和N －1次复数加法。算出全部N 点)(k X 共需2 N 次复数乘法和)1(-N N 次复数加法。即计算量是与2 N 成正比的。 FFT 的基本思想：将大点数的DFT 分解为若干个小点数DFT 的组合，从而减少运算量。 N W 因子具有以下两个特性，可使DFT 运算量尽量分解为小点数的DFT 运算： （1） 周期性：k N n N kn N n N k N W W W )()(++== （2） 对称性：k N N k N W W -=+) 2/( 利用这两个性质，可以使DFT 运算中有些项合并，以减少乘法次数。例子：求当N ＝4时，X(2)的值 ) ()]3()1([)]2()0([)()]3()1([)]2()0([)3()2()1()0()()2(0 424046 44424043 24对称性－＝周期性W x x x x W x x W x x W x W x W x W x W n x X n n +++++=+++==∑= 通过合并，使乘法次数由4次减少到1次，运算量减少。 FFT 的算法形式有很多种，但基本上可以分为两大类：按时间抽取（DIT ）和按频率抽取（DIF ）。 4.1 按时间抽取（DIT ）的FTT 为了将大点数的DFT 分解为小点数的DFT 运算，要求序列的长度N 为复合数，最常用的是M N 2=的情况（M 为正整数）。该情况下的变换称为基2FFT 。下面讨论基2情况的算法。 先将序列x(n)按奇偶项分解为两组 ?? ?=+=) ()12()()2(21r x r x r x r x 12,,1,0-=N r Λ 将DFT 运算也相应分为两组 = =)]([)(n x DFT k X ∑-=1 )(N n kn N W n x ∑∑-=-== 1 n 0 1 0)()(N n kn N N n n kn N W n x W n x 为奇数为偶数＋ ∑∑-=+-=++ 1 2/0 )12(1 2/0 2)12()2(N r k r N N r rk N W r x W r x ＝ FFT算法的DSP实现 对于离散傅里叶变换(DFT)的数字计算，FFT是一种有效的方法。一般假定输入序列是复数。当实际输入是实数时，利用对称性质可以使计算DFT非常有效。 一个优化的实数FFT算法是一个组合以后的算法。原始的2N个点的实输入序列组合成一个N点的复序列，之后对复序列进行N点的FFT运算，最后再由N点的复数输出拆散成2N点的复数序列，这2 N点的复数序列与原始的2N点的实数输入序列的DFT输出一致。 使用这种方法，在组合输入和拆散输出的操作中，FFT运算量减半。这样利用实数FFT 算法来计算实输入序列的DFT的速度几乎是一般FFT算法的两倍。下面用这种方法来实现一个256点实数FFT（2N=256）运算。 1.实数FFT运算序列的存储分配 如何利用有限的DSP系统资源，合理的安排好算法使用的存储器是一个比较重要的问题。本文中，程序代码安排在0x3000开始的存储器中，其中0x3000~0x3080存放中断向量表，FFT程序使用的正弦表和余弦表数据(.data段)安排在0xc00开始的地方，变量(.bss段定义) 存放在0x80开始的地址中。另外，本文中256点实数FFT程序的数据缓冲位0x2300~0x23ff , FFT变换后功率谱的计算结果存放在0x2200~0x22ff中。连续定位.cmd文件程序如下： MEMORY { PAGE 0: IPROG: origin = 0x3080,len=0x1F80 VECT: lorigin=0x3000,len=0x80 EPROG: origin=0x38000,len=0x8000 PAGE 1: USERREGS: origin=0x60,len=0x1c BIOSREGS: origin=0x7c,len=0x4 IDA TA: origin=0x80,len=0xB80 EDA TA: origin=0xC00,len=0x1400 } SECTIONS { .vectors: { } > VECT PAGE 0 .sysregs: { } > BIOSREGS PAGE 1 .trcinit: { } > IPROG PAGE 0 .gblinit: { } > IPROG PAGE 0 .bios: { } > IPROG PAGE 0 frt: { } > IPROG PAGE 0 .text: { } > IPROG PAGE 0 .cinit: { } > IPROG PAGE 0 .pinit: { } > IPROG PAGE 0 .sysinit: { } > IPROG PAGE 0 .data: { } > EDATA PAGE 1 .bss: { } > IDATA PAGE 1 .far: { } > IDATA PAGE 1 .const: { } > IDATA PAGE 1 .switch: { } > IDATA PAGE 1 数字信号处理期末大作业FFT的发展史、现状及典型算法 班级 学号： 姓名： FFT的发展史、现状及典型算法 傅里叶分析已有200多年的历史，目前FFT及其校正算法在工程实际中仍在广泛应用，展现了其不竭的生命力。本次作业我们论述FFT的现状，发展史以及一些算法，去详细了解、扩展这一算法，巩固所学知识。 一．FFT的简介 傅里叶变换是一种将信号从时域变换到频域的变换形式，然而当N很大的时候，求一个N点的DFT要完成N*N次复数乘法和N*（N-1）次复数加法，计算量非常大，所以人们开始探索一种简便的算法对于一个较大的N进行傅里叶变换。在20世纪60年代由Cooley和Tukey提出了快速傅里叶变换算法，它是快速计算DFT的一种简单高效的方法。 关于何为FFT，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。举个例子，设x(n)为N项的复数序列，由DFT 变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N^2次运算。当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+（N^2）/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。 所以使用FFT算法，可以大大提高傅里叶变换的运算速度，运算时间缩短一到两个数量级，从而使DFT变换应用迅速普及，不仅在频谱分析，而且在线性卷积、线性相关等方面得到广泛应用。 《DSP原理及其应用》实验设计报告 实验题目：FFT算法及滤波器的设计 摘要 随着信息科学的迅猛发展，数据采集与处理是计算机应用的一门关键技术，它主要研究信息数据的采集、存储和处理。而数字信号处理器(DSP)芯片的出现为实现数字信号处理算法提供了可能。数字信号处理器（DSP）以其特有的硬件体系结构和指令体系成为快速精确实现数字信号处理的首选工具。DSP芯片采用了哈佛结构，以其强大的数据处理功能在通信和信号处理等领域得到了广泛应用，并成为研究的热点。 本文主要研究基于TI的DSP芯片TMS320c54x的FFT算法、FIR滤波器和IIR滤波器的实现。首先大概介绍了DSP和TMS320c54x的结构和特点并详细分析了本系统的FFT变换和滤波器的实现方法。 关键词：DSP、TMS320c54x、FFT、FIR、IIR Abstract With the rapid development of information science, data acquisition and processing is a key technology of computer applications, the main research of it is collection, storage and processing of information data. The emergence of the digital signal processor (DSP) chip offers the potential for the realization of the digital signal processing algorithm. Digital signal processor (DSP), with its unique hardware system structure and instruction system become the first tool of quickly and accurately realize the digital signal processing.DSP chip adopted harvard structure, with its powerful data processing functions in the communication and signal processing, and other fields has been widely applied, and become the research hot spot. This paper mainly studies the FFT algorithm based on TMS320c54x DSP chip of TI, the realization of FIR filter and IIR filter. First introduced the DSP and TMS320c54x briefly, then analyzed in detail the structure and characteristics of the system of the realization of FFT transform and filter method. Keyword: DSP、TMS320c54x、FFT、FIR、IIR A1=str2double(get(handles.edit8,'String')); A2=str2double(get(handles.edit9,'String')); F1=str2double(get(handles.edit10,'String')); F2=str2double(get(handles.edit11,'String')); Fs=str2double(get(handles.edit12,'String')); N=str2double(get(handles.edit13,'String')); t=[0:1/Fs:(N-1)/Fs]; x=A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t); %信号x的离散值 axes(handles.axes1) %在axes1中作原始信号图 plot(x); grid on m=nextpow2(x);N=2^m; % 求x的长度对应的2的最低幂次m if length(x) 实数FFT算法的设计及其C语言实现 目前国内有关数字信号处理的教材在讲解快速傅里叶变换(FFT)时，都是以复数FFT为重点，实数FFT算法都是一笔带过，书中给出的具体实现程序多为BASIC或FORTRAN程序并且多数不能真正运行。鉴于目前在许多嵌入式系统中要用到FFT运算，如以DSP为核心的交流采样系统、频谱分析、相关分析等。本人结合自己的实际开发经验，研究了实数的FFT算法并给出具体的C语言函数，读者可以直接应用于自己的系统中。 首先分析实数FFT算法的推导过程，然后给出一种具体实现FFT算法的C语言程序，可以直接应用于需要FFT运算的单片机或DSP等嵌入式系统中。 1 倒位序算法分析 按时间抽取(DIT)的FFT算法通常将原始数据倒位序存储，最后按正常顺序输出结果X(0)，X(1)，...，X(k)，...。假设一开始，数据在数组float dataR[128]中，我们将下标i表示为(b6b5b4b3b2b1b0)b，倒位序存放就是将原来第i个位置的元素存放到第(b0b1b2b3b4b5b6)b 的位置上去.由于C语言的位操作能力很强，可以分别提取出b6、b5、b4、b3、b2、b1、b0，再重新组合成b0、b1、b2、b3、b4、b5、b6，即是倒位序的位置。程序段如下(假设128点FFT)： /* i为原始存放位置，最后得invert_pos为倒位序存放位置*/ int b0=b1=b2=b3=b4=b5=6=0; b0=i0x01; b1=(i/2)0x01; b2=(i/4)0x01; b3=(i/8)0x01; b4=(i/16)0x01; b5=(i/32)0x01; b6=(i/64)0x01; /*以上语句提取各比特的0、1值*/ invert_pos=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6; 大家可以对比教科书上的倒位序程序，会发现这种算法充分利用了C语言的位操作能力，非常容易理解而且位操作的速度很快。 2 实数蝶形运算算法的推导 我们首先看一下图1所示的蝶形图。 ARM或单片机用的FFT算法，用于信号处理。经过愚昧的本人优化，提高了计算效率 （在aduc7026 @41MHz FFT256点环境下运算时间为0.06s左右） PS：以下有两部分（fft.h和fft.c） 【复制以下内容改名为fft.h】 #ifndef __FFT_H__ #define __FFT_H__ #include 按时间抽取的基2FFT 算法分析及MATLAB 实现 一、DIT-FFT 算法的基本原理 基2FFT 算法的基本思想是把原始的N 点序列依次分解成一系列短序列，充分利用旋转因子的周期性和对称性，分别求出这些短序列对应的DFT ，再进行适当的组合，得到原N 点序列的DFT ，最终达到减少运算次数，提高运算速度的目的。 按时间抽取的基2FFT 算法，先是将N 点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n)，再分别进行DFT 运算，求出与之对应的X1(k)和X2(k)，然后利用图1所示的运算流程进行蝶形运算，得到原N 点序列的DFT 。只要N 是2的整数次幂，这种分解就可一直进行下去，直到其DFT 就是本身的1点时域序列。 图1 DIT-FFT 蝶形运算流图 二、DIT-FFT 算法的运算规律及编程思想 1.原位计算 对N=M 2点的FFT 共进行M 级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，经过M 级运算后，原来存放输入序列数据的N 个存储单元中可依次存放X(k)的N 个值，这种原位(址)计算的方法可节省大量内存。 2.旋转因子的变化规律 N 点DIT ―FFT 运算流图中，每个蝶形都要乘以旋转因子p W N ，p 称为旋转因子的指数。例如N ＝8 ＝3 2 时各级的旋转因子： 第一级：L=1， 有1个旋转因子：p W N =J /4W N =J 2L W J=0 第二级：L=2，有2个旋转因子：p W N =J /2W N =J 2L W J=0,1 第三级：L=3，有4个旋转因子：p W N =J W N =J 2L W J=0,1,2,3 对于N ＝M 2的一般情况，第L 级共有1 -L 2 个不同的旋转因子： p W N =J 2L W J=0,1,2,… ,1 -L 2－1 L 2=M 2×M -L 2 = N ·M -L 2 故： 按照上面两式可以确定第L 级运算的旋转因子 数字通信课程设计报告书 课题名称 基于matlab 的FFT 算法程序设计 姓 名 学 号 院 系 物理与电信工程系 专 业 电子信息工程 指导教师 2010年 01 月15日 ※※※※※※※※※ ※ ※ ※※ ※※ ※※ ※※※※※ ※※ 2007级数字通信 课程设计 基于matlab的FFT算法程序设计 0712401-36 李晔 （湖南城市学院物理与电信工程系通信工程专业,益阳,413000） 一、设计目的 1.通过该设计，进一步了解MATLAB软件。 2.通过该设计，进一步熟悉MATLAB的语法规则和编辑方式。 3.通过该设计，掌握傅里叶变换的含义和方法。 二、设计的主要要求 掌握Fourier变换，解了关于MATLAB软件在数字信号处理方面的应用，熟悉MATLAB的语法规则和编程。用MATLAB实现快速Fourier变换。 三、整体设计方案 对信号x=sin(2*pi*f0*t)进行频谱分析，用MATLAB仿真。选取抽样频率为fs=100Hz,依照下列条件用MATLAB软件对信号xt进行傅里叶变换y=fft(xt,N)并绘制频谱图，观察所产生的六幅频谱图进行对比，并进行分析。 四、程序设计 fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率 %生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(321); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 m=length(y); f=(0:m/2-1)'*fs/m;%进行对应的频率转换 figure(1); subplot(322); plot(f,mag(1:m/2));%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值'); title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱 sq=abs(y); figure(1); subplot(323); plot(f,sq(1:m/2)); xlabel('频率(Hz)'); ylabel('均方根谱'); 毕业设计 [论文] 题目：FFT算法研究及基2-FFT算法的C语言实现学院：电气与信息工程学院 专业：电气工程及其自动化 姓名：XXX 学号：XXXXXX 指导老师：XXX 完成时间：2015年06月01日 摘要 离散傅立叶变换（DFT）常常用于计算信号处理。DFT算法可以得到信号的频域特性，因为该算法在计算上是密集的，长时间的使用时，计算机不能实时进行信号处理。所以DFT被发现之后的相当长时间内是没被应用到实际的项目。到了二十世纪六十年代中期一种新的计算方法被研究者发现出来，它就是FFT。FFT 并不是一种新的获取频域特征的方式，而是离散傅里叶变换的一种快速实现算法。 数字信号处理在当今科技发展中发展很迅速，不但是在传统的通信领域，其他方面也经常用到。利用快速傅里叶变换，实现了信号频域的变换处理。对于信号的处理，往往会和数学中的算法联系到一起。如果处理得当，将会对气象，地理信息等的发展，起到举足轻重的作用，同时对世界其他领域的发展有很大的促进作用。 关键词: FFT算法，C语言，编译实现 Abstract Discrete Fourier Transform (DFT) is often used to calculate the signal processing to obtain frequency domain signals. DFT algorithm can get the frequency domain characteristics of the signal, because the algorithm is computationally intensive, long-time use, the computer is not conducive to real-time signal processing. So DFT since it was discovered in a relatively long period of time is not to be applied to the actual projects until a new fast discrete Fourier calculation method --FFT is found in discrete Fourier transform was able to actually project has been widely used. FFT is not a new way to get the frequency domain, but the discrete Fourier transform of a fast algorithm. Fast Fourier Transform (FFT) is a digital signal processing important tool that the time domain signal into a frequency-domain signal processing. matched filtering has important applications. FFT is a discrete Fourier transform (DFT) is a fast algorithm, it can be a signal from the time domain to the frequency domain. Some signals in the time domain is not easy to observe the characteristics of what is, but then if you change the signal from the time domain to the frequency domain, it is easy to see out of. This design is required to be familiar with the basic principles of FFT algorithm, based on the preparation of C language program to achieve FFT algorithm real number sequence. Keywords: FFT algorithm, C language compiler to achieve 基2的DIT蝶形算法源代码及注释如下： /************FFT***********/ //整个程序输入和输出利用同一个空间x[N]，节约空间 #include 实验二快速傅立叶变换（FFT）算法实验 一．实验目的 1．加深对DFT算法原理和基本性质的理解； 2．熟悉FFT算法原理和FFT子程序的应用； 3．学习用FFT对连续信号和时域信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。 二．实验设备 计算机，CCS 2.0 版软件，实验箱，DSP仿真器，短接块，导线。 三．基本原理 1．离散傅立叶变换DFT的定义：将时域的采样变换成频域的周期性离散函数，频域的采样也可以变换成时域的周期性离散函数，这样的变换称为离散傅立叶变换，简称DFT。 2．FFT是DFT的一种快速算法，将DFT的N2步运算减少为（N/2）log2N步，极大的提高了运算的速度。 3．旋转因子的变化规律。 4．蝶形运算规律。 5．基2FFT算法。 四．实验步骤 1．复习DFT的定义、性质和用DFT作谱分析的有关内容； 2．复习FFT算法原理与编程思想，并对照DIT-FFT运算流程图和程序框图，了解本实验提供的FFT子程序； 3．阅读本实验所提供的样例子程序； 4．运行CCS软件，对样例程序进行跟踪，分析结果；记录必要的参数。 5．填写实验报告。 6．提供样例程序实验操作说明 1）实验前的准备 “语音处理单元”的拨码开关设置： 在信号源单元中，设置左路信号源产生低频正弦波信号，右路产生高频正弦波信号。实验箱上电，用示波器分别观测OUT1和OUT2输出的模拟信号，并调节电位器直至低频正弦波信号为100Hz/1V左右；高频正弦波信号为6KHz/1V左右；将S3中的拨码开关2打到ON，用示波器观测OUT1输出的混叠信号波形。 用导线连接“信号源单元”中2号孔接口OUT1和语音处理单元中的2号孔接口“IN”； 正确完成计算机、DSP仿真器和实验箱的连接后，系统上电. 2）实验过程 启动CCS 2.0，用Project/Open打开“ExpFFT01.pjt”工程文件；双击“ExpFFT01.pjt”及“Source”可查看各源程序；加载“ExpFFT01.out”； 在主程序中，k++处设置断点；单击“Run”运行程序，或按F5运行程序；程序将运行 利用FFT 计算卷积 一．线卷积的作用及定义 线卷积包括卷积积分和卷积和。 1．线卷积的作用 求解线性系统对任意激励信号的零态响应。 2．卷积积分 ) (*)(d )()()(t h t x t h x t y =-= ? ∞∞ -τττ 3．卷积和 离散系统的时域分析是，已知离散系统的初始状态和输入信号（激励），求离散系统的输出（响应），两种方法：递推解法和离散卷积法。 卷积和：)()()()()(n h n x m n h m x n y m *=-= ∑ ∞ -∞ = 二．圆周卷积的定义 圆周移位：一周期为N 的周期序列, 可视为一主值序列在圆周上的循环移位。周期序列在时间轴上左移 右移m 反时针 转称为圆周移位。 时域圆周卷积（循环卷积） )()()(n h n x n y ?=()()()∑ -=-= 1 )(N m N N n R m n h m x 条件：两序列实现圆卷积的条件是：长度相等，如果不相等, 可通过增补零值来使之相等。 特点：卷积求和范围只在10-≤≤N m 有限区间进行；卷积时不作反褶平移, 而是反褶圆移 步骤：量置换→反褶→圆移→相乘→求和。 三．两者的关系 有限长序列的圆卷积和线卷积的关系 在一般情况下，两序列的圆卷积和线卷积是不相等的，这是因为：线卷积是 平移, 结果长度为121-+=N N L ；而圆卷积是圆移，结果长度为2 1 N N L ==。只有 在两卷积的结果长度相时，二者才有相同的结果。解决方法是：在作圆卷积时，通过加零的方法，使两序列的长度都增加到121-+=N N L ，此时，圆卷积的结果和线卷积同。 四．利用FFT 计算卷积 工程实际需要解决的卷积：)()()(n h n x n y *=，但其计算量很大。 而圆卷积为：)()()(n h n x n y ?=，便于采用FFT 算法, 故计算速度快。若将线卷积的两个序列用增补零的方法将长度取为一致，此时两序列的离散线卷积和圆周卷积结果是相等的，这样就则可以通过圆卷积来快速计算线卷积。 1、 利用FFT 计算卷积的步骤 （1）设两序列原长度分别为：N 和M ，将长度增加到1-+≥M N L （L 为2的整数次幂）； （2）用FFT 法求加长序列的DFT 频谱； （3）计算两序列DFT 频谱的乘积； （4）用IFFT 求DFT 频谱乘积的逆变换，便得两序列的离散线卷积。 2、分段快速卷积 设)(n x 为长序列，)(n h 为短序列，长度为M ，则两序列的离散线卷积可以写成如 下 形 式 ， ∑∑∑-=-+=-=+-+ +-+ -= *=1 1 )1(1 2)()()()()()()()()(N m n K kN m N N m m N h m x m N h m x m N h m x n h n x n y 上述每个子段长度为N 。为便于圆卷积计算，将长度通过补零加长为：1-+=M N L x (n 0 n h (n 根据各子段()n x k 增补零的部位不一样而分两种算法。 实验2 基2时域抽选的FFT 程序设计与调试 一、实验目的 掌握信号处理，尤其是数字信号处理的基本原理和方法。要求能通过实验熟练掌握基2时域抽选的快速傅立叶变换算法（FFT ）的基本原理，了解二维及多维快速傅立叶变换算法。 二、实验原理 1．复数类型 对于FFT 算法涉及的复数运算，使用自定义的COMPLEX 来定义复数类型，其使用方法与常规类型（如int,float,double ）相似。 typedef struct { float real, imag; } COMPLEX; 2．FFT 基本原理 FFT 改进了DFT 的算法，减少了运算量，主要是利用了旋转因子W 的两个性质： （a ）W 的周期性：W = W （b) W 的对称性：W =-W FFT 把N 点DFT 运算分解为两组N/2点的DFT 运算，然后求和： )()()(21k X W k X k X k N += 1,,1,0 ),()()2 (2 21-=-=+ N k N k k X W k X N k X 其中， ∑∑∑∑-=-=-=-=+== = = 1 1 2 21 1 112 2 2 2 2 2 2 2 )12()()()2()()(N N N N N N N N r rk r rk r rk r rk W r x W r x k X W r x W r x k X 在计算X 1(k)与X 2(k)时，仍利用上述公式，把它们看成是新的X(k)。如此递归下去，便是FFT 算法。 3．蝶形运算 从基2时域抽选FFT 运算流图可知： ① 蝶形两节点的距离为2m-1，其中，m 表示第m 列，且m =1,… ,L 。 例如N=8=23， 第一级(列)距离为21-1=1， 第二级(列)距离为22-1=2， 第三级(列)距离为23-1=4。 ② 考虑蝶形运算两节点的距离为2m-1，蝶形运算可表为： X m (k)=X m-1(k)+X m-1(k+2m-1) W N r X m (k+2m-1)= X m-1(k)-X m-1(k+2m-1) W N r 由于N 为已知，所以将r 的值确定即可确定W N r 。为此，令k=(n 2n 1n 0)2 ,再将k 左移(L-m)位，右边位置补零，就可得到(r)2 的值，即(r)2 =(k)22L-m 。 例如 N=8=23FFT 算法详解
FFT算法的意义
按时间抽取的基2FFT算法分析
FFT算法的DSP实现
fft的发展、现状、典型算法
FFT算法及IIR、FIR滤波器的设计
用matlab实现fft算法
实数FFT算法的设计及其C语言实现
FFT算法(查表法)
按时间抽取的基2FFT算法分析及MATLAB实现
基于matlab的FFT算法程序设计
FFT算法研究及基2-FFT算法的C语言实现
fft算法代码注释及流程框图
快速傅立叶变换(FFT)算法实验
利用FFT计算卷积
信号处理 FFT算法