FFT算法的意义
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。如果要提高频率
分辨力，则必须增加采样点数，也即采样时间。频率分辨率和
采样时间是倒数关系。
假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是
An=根号a*a+b*b，相位就是Pn=atan2(b,a)。根据以上的结果，
就可以计算出n点（n≠1，且n<=N/2）对应的信号的表达式为：
An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。
对于n=1点的信号，是直流分量，幅度即为A1/N。
由于FFT结果的对称性，通常我们只使用前半部分的结果，
即小于采样频率一半的结果。
好了，说了半天，看着公式也晕，下面圈圈以一个实际的
信号来做说明。
假设我们有一个信号，它含有2V的直流分量，频率为50Hz、
相位为-30度、幅度为3V的交流信号，以及一个频率为75Hz、
相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下：
S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)
式中cos参数为弧度，所以-30度和90度要分别换算成弧度。
我们以256Hz的采样率对这个信号进行采样，总共采样256点。
按照我们上面的分析，Fn=(n-1)*Fs/N，我们可以知道，每两个
点之间的间距就是1Hz，第n个点的频率就是n-1。我们的信号
有3个频率：0Hz、50Hz、75Hz，应该分别在第1个点、第51个点、第76个点上出现峰值，其它各点应该接近0。实际情况如何呢？
我们来看看FFT的结果的模值如图所示。
图1 FFT结果
从图中我们可以看到，在第1点、第51点、和第76点附近有比较大的值。我们分别将这三个点附近的数据拿上来细看：
1点：512+0i
2点：-2.6195E-14 - 1.4162E-13i
3点：-2.8586E-14 - 1.1898E-13i
50点：-6.2076E-13 - 2.1713E-12i
51点：332.55 - 192i
52点：-1.6707E-12 - 1.5241E-12i
75点：-2.2199E-13 -1.0076E-12i
76点：3.4315E-12 + 192i
77点：-3.0263E-14 +7.5609E-13i
很明显，1点、51点、76点的值都比较大，它附近的点值
都很小，可以认为是0，即在那些频率点上的信号幅度为0。
接着，我们来计算各点的幅度值。分别计算这三个点的模值，
结果如下：
1点：512
51点：384
76点：192
按照公式，可以计算出直流分量为：512/N=512/256=2；
50Hz信号的幅度为：384/(N/2)=384/(256/2)=3；75Hz信号的
幅度为192/(N/2)=192/(256/2)=1.5。可见，从频谱分析出来
的幅度是正确的。
然后再来计算相位信息。直流信号没有相位可言，不用管
它。先计算50Hz信号的相位，atan2(-192, 332.55)=-0.5236,
结果是弧度，换算为角度就是180*(-0.5236)/pi=-30.0001。再
计算75Hz信号的相位，atan2(192, 3.4315E-12)=1.5708弧度，
换算成角度就是180*1.5708/pi=90.0002。可见，相位也是对的。
根据FFT结果以及上面的分析计算，我们就可以写出信号的表达
式了，它就是我们开始提供的信号。
总结：假设采样频率为Fs，采样点数为N，做FFT之后，某
一点n（n从1开始）表示的频率为：Fn=(n-1)*Fs/N；该点的模值除以N/2就是对应该频率下的信号的幅度（对于直流信号是除以N）；该点的相位即是对应该频率下的信号的相位。相位的计算
可用函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角
度值，范围从-pi到pi。要精确到xHz，则需要采样长度为1/x秒
的信号，并做FFT。要提高频率分辨率，就需要增加采样点数，
这在一些实际的应用中是不现实的，需要在较短的时间内完成
分析。解决这个问题的方法有频率细分法，比较简单的方法是
采样比较短时间的信号，然后在后面补充一定数量的0，使其长度
达到需要的点数，再做FFT，这在一定程度上能够提高频率分辨力。具体的频率细分法可参考相关文献。
[附录：本测试数据使用的matlab程序]
close all; %先关闭所有图片
Adc=2; %直流分量幅度
A1=3; %频率F1信号的幅度
A2=1.5; %频率F2信号的幅度
F1=50; %信号1频率(Hz)
F2=75; %信号2频率(Hz)
Fs=256; %采样频率(Hz)
P1=-30; %信号1相位(度)
P2=90; %信号相位(度)
N=256; %采样点数
t=[0:1/Fs:N/Fs]; %采样时刻
%信号
S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180); %显示原始信号
plot(S);
title('原始信号');
figure;
Y = fft(S,N); %做FFT变换
Ayy = (abs(Y)); %取模
plot(Ayy(1:N)); %显示原始的FFT模值结果
title('FFT 模值');
figure;
Ayy=Ayy/(N/2); %换算成实际的幅度
Ayy(1)=Ayy(1)/2;
F=([1:N]-1)*Fs/N; %换算成实际的频率值
plot(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果
title('幅度-频率曲线图');
figure;
Pyy=[1:N/2];
for i=1:N/2
Pyy(i)=phase(Y(i)); %计算相位
Pyy(i)=Pyy(i)*180/pi; %换算为角度
end;
plot(F(1:N/2),Pyy(1:N/2)); %显示相位图
title('相位-频率曲线图');
FFT 算法详解
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。 现在圈圈就根据实际经验来说说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。如果要提高频率
实验二 FFT算法的MATLAB实现
班级：学号：姓名 实验二FFT算法的MATLAB实现 （一）实验目的： （1）掌握用matlab进行FFT在数字信号处理中的高效率应用。 （2）学习用FFT对连续信号和时域离散信号进行谱分析。 （二）实验内容及运行结果： 题1：若x(n)=cos(nπ/6)是一个N=12的有限序列，利用MATLAB计算它的DFT 并进行IDFT变换同时将原图与IDFT变换后的图形进行对比。当求解IFFT变换中，采样点数少于12时，会产生什么问题。 程序代码： N=12; n=0:11; Xn=cos(n*pi/6); k=0:11; nk=n'*k; WN=exp(-j*2*pi/N) WNnk=WN.^nk XK=Xn*WNnk; figure(1) stem(Xn) figure(2) stem(abs(XK)) 运行结果：
IFFT变换中，当采样点数少于12时图像如下图显示：
分析：由图像可以看出，当采样点数小于12时，x(n)的频谱不变，周期为6，而XK 的频谱图发生改变。 题2：对以下序列进行谱分析 132()()103()8470x n R n n n x n n n =+≤≤?? =-≤≤??? 其他n 选择FFT 的变换区间N 为8和16点两种情况进行频谱分析，分别打印其幅频特 性曲线并进行对比、分析和讨论。 ㈠ 程序代码： x=ones(1,3);nx=0:2; x1k8=fft(x,8); F=(0:length(x1k8)-1)'*2/length(x1k8); %进行对应的频率转换 stem(f,abs(x1k8));%8点FFT title('8点FFTx_1(n)'); xlabel('w/pi'); ylabel('幅度'); N=8时：
按时间抽取的基2FFT算法分析及MATLAB实现
按时间抽取的基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 级运算的旋转因子 \ 3、同一级中，同一旋转因子对应蝶形数目 第L 级FFT 运算中，同一旋转因子用在L -M 2 个蝶形中； 4、同一级中，蝶形运算使用相同旋转因子之间相隔的“距离” 第L 级中，蝶距：D=L 2； 5、同一蝶形运算两输入数据的距离 在输入倒序，输出原序的FFT 变换中，第L 级的每一个蝶形的2个输入数据相距：B=1 -L 2。 6、码位颠倒 输入序列x(n)经过M 级时域奇、偶抽选后，输出序列X(k)的顺序和输入序列的顺序关系为倒位关系。 '
按时间抽取的基2FFT算法分析
第四章 快速傅里叶变换 有限长序列可以通过离散傅里叶变换(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实现
1 引言 CCS(Code Composer Studio)是TI公司的DSP集成开发环境。它提供了环境配置、源文件编辑、程序调试、跟踪和分析等工具,帮助用户在一个软件环境下完成编辑、编译链接、调试和数据分析等工作。与TI提供的早期软件开发工具相比,利用CCS能够加快软件开发进程,提高工作效率。CCS一般工作在两种模式下:软件仿真器和与硬件开发板相结合的在线编程。前者可以脱离DSP芯片,在PC机上模拟DSP指令集与工作机制,主要用于前期算法实现和调试。后者实时运行在DSP芯片上,可以在线编制和调试应用程序。 2 C语言和汇编语言的混合编程 TMS320 C5000系列的软件设计通常有三种方法: (1) 用C语言开发; (2) 用汇编语言开发; (3) C和汇编的混合开发。 其中用C语言开发具有兼容性和可移植的优点,有利于缩短开发周期和减少开发难度,但是在运算量较大的情况下,C代码的效率还是无法和手工编写的汇编代码的效率相比,比如FFT运算,用汇编语言开发的效率高,程序执行速度快,而且可以合理利用芯片的硬件资源,但是开发难度较大,开发周期长,而且可读性和可移植性差。C和汇编的混合编程则可以充分利用前两者的优点,以达到最佳利用DSP资源的目的。但是,采用C和汇编语言混合编程必须遵循相关函数调用规则和寄存器调用规则,否则会给程序的开发带来意想不到的问题。 2.1 C语言和汇编语言混合编程的四种方法 (1) 独立编写汇编程序和C程序,分开编译或汇编成各自的目标代码模块,再用链接器将二者链接起来。这种方法比较灵活,但是设计者必须自己维护各汇编模块的入口和出口代码,自己计算传递的参数在堆栈中的偏移量,工作量较大,但是能做到对程序的绝对控制。 (2) 在C程序中使用汇编程序中定义的变量和常数。 (3) 在C程序中内嵌汇编语句。这种方法可以实现C语言无法实现的一些硬件控制功能,如修改中断控制寄存器。 (4) 将C语言编译生成相应的汇编代码,手工修改和优化C编译器生成的汇编代码。采用这种方法可以控制C编译器,从而产生具有交叉列表的汇编程序,而设计者可以对其中的汇编语句进行修改,然后对汇编程序进行编译,产生目标文件。
FFT-C快速傅里叶变换超级详细的原代码
快速傅立叶变换（FFT）的C++实现收藏 标准的离散傅立叶DFT 变换形式如： y k=Σj=0n-1a jωn-kj = A (ωn-k). （ωn k为复数1 的第k 个n 次方根，且定义多项式A (x)=Σj=0n-1a j x j） 而离散傅立叶逆变换IDFT (Inverse DFT)形式如：a j=(Σk=0n-1y kωn kj)/n . yk=Σj=0n-1 ajωn-kj = A (ωn-k). （ωnk 为复数1 的第k 个n 次方根，且定义多项式 A (x) = Σj=0n-1 ajxj ） 而离散傅立叶逆变换IDFT (Inverse DFT)形式如：aj=(Σk=0n-1 ykωnkj)/n . 以下不同颜色内容为引用并加以修正： 快速傅立叶变换（Fast Fourier Transform，FFT）是离散傅立叶变换（Discrete Fourier transform，DFT）的快速算法，它是根据离散傅立叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。 设Xn 为N 项的复数序列，由DFT 变换，任一Xi 的计算都需要N 次复数乘法和N -1 次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N 项复数序列的Xi ，即N 点DFT 变换大约就需要N2 次运算。当N =1024 点甚至更多的时候，需要N2 = 1048576 次运算，在FFT 中，利用ωn 的周期性和对称性，把一个N 项序列（设N 为偶数），分为两个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 变换就只需要N * log2N 次的运算，N = 1024 点时，运算量仅有10240 次，是先前的直接算法的1% ，点数越多，运算量的节约就越大，这就是FFT 的优越性。 FFT 的实现可以自顶而下，采用递归，但是对于硬件实现成本高，对于软件实现都不够高效，改用迭代较好，自底而上地解决问题。感觉和归并排序的迭代版很类似，不过先要采用“位反转置换”的方法把Xi 放到合适的位置，设i 和j 互为s = log2N 位二进制的回文数，假设s = 3, i = (110)2 = 6, j = (011)2 = 3, 如果i ≠j , 那么Xi 和Xj 应该互换位置。（关于这个回文数的生成，是很有趣而且是很基本的操作，想当初偶初学C++ 的时候就有这样的习题。）当“位反转置换”完成后，先将每一个Xi 看作是独立的多项式，然后两个两个地将它们合并成一个多项式（每个多项式有2 项），合并实际上是“蝶形运算”（Butterfly Operation, 参考《算法导论》吧^_^），继续合并（第二次的每个多项式有4 项），直到只剩
fft的发展、现状、典型算法
数字信号处理期末大作业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变换应用迅速普及，不仅在频谱分析，而且在线性卷积、线性相关等方面得到广泛应用。
用matlab实现fft算法
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) 《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 实数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 实验2 FFT 算法实现 2.1 实验目的 1、 加深对快速傅里叶变换的理解。 2、 掌握FFT 算法及其程序的编写。 3、 掌握算法性能评测的方法。 2.2 实验原理 一个连续信号)(t x a 的频谱可以用它的傅立叶变换表示为 dt e t x j X t j a a Ω-+∞ ∞-?= Ω)()( （2-1） 如果对该信号进行理想采样，可以得到采样序列 )()(nT x n x a = （2-2） 同样可以对该序列进行z 变换，其中T 为采样周期 ∑+∞∞--=n z n x z X )()( （2-3） 当ωj e z =的时候，我们就得到了序列的傅立叶变换 ∑+∞∞-=n j j e n x e X ωω)()( （2-4） 其中ω称为数字频率，它和模拟域频率的关系为 s f T /Ω=Ω=ω （2-5） 式中的s f 是采样频率。上式说明数字频率是模拟频率对采样率s f 的归一化。同模拟域的情况相似，数字频率代表了序列值变化的速率，而序列的傅立叶变换称为序列的频谱。序列的傅立叶变换和对应的采样信号频谱具有下式的对应关系。 ∑+∞∞--=)2(1)(T m j X T e X a j πωω （2-6） 即序列的频谱是采样信号频谱的周期延拓。从式（2-6）可以看出，只要分析采样序列的频谱，就可以得到相应的连续信号的频谱。注意：这里的信号必须是带限信号，采样也必须满 足Nyquist 定理。 在各种信号序列中，有限长序列在数字信号处理中占有很重要的地位。无限长的序列也往往可以用有限长序列来逼近。对于有限长的序列我们可以使用离散傅立叶变换（DFT ），这一变换可以很好地反应序列的频域特性，并且容易利用快速算法在计算机上实现当序列的长度是N 时，我们定义离散傅立叶变换为： ∑-===10)()]([)(N n kn N W n x n x DFT k X （2-7） 其中N j N e W π 2-=，它的反变换定义为： ∑-=-==10)(1)]([)(N k kn N W k X N k X IDFT n x （2-8） 根据式（2-3）和（2-7）令k N W z -=，则有 ∑-====-10)]([)(|)(N n nk N W z n x DFT W n x z X k N （2-9） 可以得到k N j k N e W z z X k X π2|)()(===-,k N W -是z 平面单位圆上幅角为k N πω2=的点，就是将单位圆进行N 等分以后第k 个点。所以，X(k)是z 变换在单位圆上的等距采样，或者说是序列傅立叶变换的等距采样。时域采样在满足Nyquist 定理时，就不会发生频谱混淆；同样地，在频率域进行采样的时候，只要采样间隔足够小，也不会发生时域序列的混淆。 DFT 是对序列傅立叶变换的等距采样，因此可以用于序列的频谱分析。在运用DFT 进行频谱分析的时候可能有三种误差，分析如下： （1）混淆现象 从式（2-6）中可以看出，序列的频谱是采样信号频谱的周期延拓，周期是2π/T ，因此当采样速率不满足Nyquist 定理，即采样频率T f s /1=小于两倍的信号（这里指的是实信号）频率时，经过采样就会发生频谱混淆。这导致采样后的信号序列频谱不能真实地反映原信号的频谱。所以，在利用DFT 分析连续信号频谱的时候，必须注意这一问题。避免混淆现象的唯一方法是保证采样的速率足够高，使频谱交叠的现象不出现。这就告诉我们，在确定信号的采样频率之前，需要对频谱的性质有所了解。在一般的情况下，为了保证高于折叠频率的分量不会出现，在采样之前，先用低通模拟滤波器对信号进行滤波。 （2）泄漏现象 实际中的信号序列往往很长，甚至是无限长序列。为了方便，我们往往用截短的序列来近似它们。这样可以使用较短的DFT 来对信号进行频谱分析。这种截短等价于给原信号序列乘以一个矩形窗函数。而矩形窗函数的频谱不是有限带宽的，从而它和原信号的频谱进行卷积以后会扩展原信号的频谱。值得一提的是，泄漏是不能和混淆完全分离开的，因为泄露导致频谱的扩展，从而造成混淆。为了减小泄漏的影响，可以选择适当的窗函数使频谱的扩散减到最小。 （3）栅栏效应 因为DFT 是对单位圆上z 变换的均匀采样，所以它不可能将频谱视为一个连续函数。 STM32F103 12-15元左右 本文将以一个实例来介绍如何使用STM32提供的DSP库函数进行FFT。 1.FFT运算效率 使用STM32官方提供的DSP库进行FFT，虽然在使用上有些不灵活（因为它是基4的FFT，所以FFT的点数必须是4^n），但其执行效率确实非常高效，看图1所示的FFT运算效率测试数据便可见一斑。该数据来自STM32 DSP库使用文档。 图1 FFT运算效率测试数据 由图1可见，在STM32F10x系列处理器上，如果使用72M的系统主频，进行64点的FFT运算，仅仅需要0.078ms而已。如果是进行1024点的FFT运算，也才需要2.138ms。 2.如何使用STM32提供的DSP库函数 2.1下载STM32的DSP库 大家可以从网上搜索下载得到STM32的DSP库，这里提供一个下载的地址：https://https://www.wendangku.net/doc/2e3980355.html,/public/STe2ecommunities/mcu/Lists/cortex_ mx_stm32/DispForm.aspx?ID=30831&RootFolder=%2fpublic%2fST e2ecommunities%2fmcu%2fLists%2fcortex%5fmx%5fstm32%2fST M32F10x%20DSP%20library%2c%20where%20is%20it 2.2添加DSP库到自己的工程项目中 下载得到STM32的DSP库之后，就可以将其添加到自己的工程项目中了。 其中，inc文件夹下的stm32_dsp.h和table_fft.h两个文件是必须添加的。stm32_dsp.h是STM32的DSP库的头文件。 src文件夹下的文件可以有选择的添加（用到那个添加那个即可）。因为我只用到了256点的FFT，所以这里我只添加了cr4_fft_256_stm32.s文件。添加完成后的项目框架如图2所示。 图2 项目框架 2.3模拟采样数据 根据采样定理，采样频率必须是被采样信号最高频率的2倍。这里，我要采集的是音频信号，音频信号的频率范围是20Hz到20KHz，所以我使用的采用频率是44800Hz。那么在进行256点FFT时，将得到44800Hz / 256 = 175Hz的频率分辨率。 为了验证FFT运算结果的正确性，这里我模拟了一组采样数据，并将该采样数据存放到了long类型的lBufInArray数组中，且该数组中每个元素的高16 位存储采样数据的实部，低16位存储采样数据的虚部(总是为0)。 为什么要这样做呢？是因为后面要调用STM32的DSP库函数，需要传入的参数规定了必须是这样的数据格式。 下面是具体的实现代码： 1 /****************************************************************** 2函数名称:InitBufInArray() 3函数功能:模拟采样数据，采样数据中包含3种频率正弦波(350Hz，8400Hz，18725Hz) 4参数说明: 5备注:在lBufInArray数组中，每个数据的高16位存储采样数据的实部， 6低16位存储采样数据的虚部(总是为0) 7作者:博客园依旧淡然（https://www.wendangku.net/doc/2e3980355.html,/menlsh/） 8 *******************************************************************/ 128点FFT 算法设计方法 X(k)=127 n 0x n =∑（）w nk 128, X （k ）=63n 0x n =∑（2）w 1282nk +63n 0x n =∑（2+1） w 128(2n+1)k =63n 0x n =∑（2）w 64nk +(63n 0x n =∑（2+1） w 64nk )w 1281k =H(k)+G(k) w 1281k , H(k)= 63 n 0h n =∑（）w 64nk ,h(n)=x(2n) G(k)= 63 n 0g n =∑（） w 64nk ,g(n)=x(2n+1) H(k) = 63n 0h n =∑（） w 64nk =7a 0=∑(7b 0h a b =∑（+8）)w 64(a+8b)k H(k)=H(k 0+8k 1)= 7a 0=∑(7b 0 h a b =∑（+8）)w 64(a+8b)(k0+8k1) =7a 0=∑(7b 0 h a b =∑（+8）w 648bk0)w 64a(k0+8k1)w 6464bk1 =7a 0 =∑(7b 0h a b =∑（+8）w 648bk0)w 64a(k0+8k1) =7a 0 =∑(7b 0h a b =∑（+8）w 8bk0)w 64a(k0+8k1), [w 6464bk1=1] =7a 0 =∑(7b 0f b =∑（）w 8bk0)w 64a(k0+8k1) =7a 0=∑(7b 0 f b =∑（）w 8bk0w 64ak0 )w 648ak1) =7a 0=∑(7 b 0f b =∑（）w 8bk0w 64ak0 )w 8ak1) w 8bk0旋转因子是1、-1、 等，可以用加减等运算实现，只有w 64ak0要乘旋转因子。 7b 0f b =∑（） w 8bk0用一个蝶8运算。 F(k0)= 7 b 0f b =∑（）w 8bk0,f(b)=h(a+2b), 红色旋转因子还未找到处理方法，使第二级也变为pe8运算 18W =134689222222 ------＝＋＋＋＋＋ 数字通信课程设计报告书 课题名称 基于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算法及IIR、FIR滤波器的设计
实数FFT算法的设计及其C语言实现
FFT算法(查表法)
实验2FFT算法实现
(完整word版)stm32F103进行FFT算法教程
128点FFT算法设计方法
基于matlab的FFT算法程序设计