基于MATLAB的FFT算法实现(论文)
基于MATLAB的FFT算法实现
摘要
MATLAB软件是目前全世界范围内非常流行的具有很强的科学计算和图形界面的软件系统。利用MATLAB的强大运算功能，可以解决数字信号处理过程中遇到的许多问题。本文给出了基于MATLAB软件实现信号DFT变换和FFT频谱分析的方法。利用MATLAB软件方法，使得设计方便、快捷，大大减轻了工作量。并且，在信号DFT变换中可以清楚得看到DFT变换结果和截取长度之间的关系。通过编程仿真可以得到序列的幅频特性曲线，便于对信号进行谱分析。
FFT（Fast Fourier Transformation），即为快速傅立叶变换，是离散傅立叶变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。
在实际应用中，FFT是最常见的数字信号处理算法，它在各种数字信号处理系统中扮演重要的角色。在信号处理过程中。频域分析往往比时域分析方便和高效，FFT是时域和频域转换的基本运算。
关键词：FFT算法；MATLAB；数字信号处理；频谱分析
THE FFT ALGORITHM BASED ON MATLAB
ABSTRACT
MATLAB software is very popular around the world have a strong scientific computing and graphic interface of the software system. Using the powerful operation function of MATLAB, can solve many problems encountered in the process of digital signal processing. In this paper, based on DFT MATLAB software to realize signal transform and FFT spectrum analysis method. Using MATLAB software method, makes the design of convenient, quick, greatly reduce the workload. And the signal DFr transform can be clearly seen in the DFT transform results and clipping of the relationship between the length. Sequences can be obtained by programming the simulation of the amplitude frequency characteristic curve, facilitate the signal spectrum analysis. FFT (Fast Fourier changed), which is Fast Fourier transform, is a Fast algorithm of discrete Fourier transform, it is according to the odd and even of discrete Fourier transform, virtual and real features, the discrete Fourier transform algorithm was improved. It the theory of Fourier transform and found nothing new, but the application in computer systems or digital system discrete Fourier transform, can be said to be into a big step.
In practical applications, the FFT is the most common form of digital signal processing algorithm, it is play an important role in all kinds of digital signal processing system. In the process of signal processing. Frequency domain analysis than time-domain analysis is convenient and efficient, FFT is the basic operation of the time domain and frequency domain transformation.
Key words:FFT algorithm; MATLAB; Digital signal processing; Spectrum analysis
目录
1 绪论 (1)
1.1 课程设计内容 (1)
1.2课程设计要求 (1)
1.3 课程设计目的 (1)
1.4 课程设计平台 (1)
2 概要设计 (3)
2.1 MATLAB及MATLAB在数字信号领域的应用…………………………………. .3
2.2 FFT算法的实现原理 (3)
3 基于MATLAB的FFT算法设计原理 (5)
3.1 总体设计流程图 (5)
3.2 语音信号的采集 (5)
3.3 语音信号的时频分析 (6)
3.4 快速傅里叶变换 (8)
3.4.1 FFT的运算规律 (10)
3.4.2 基于MATLAB的FFT所编写程序的框图 (14)
4 调试与结果 (15)
5 设计总结与心得体会 (17)
参考文献 (19)
附录 (20)
1 绪论
1.1课程设计内容
录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；在MATLAB环境下编写基2 DIT-FFT算法；利用自己编写的算法对已采集的语音信号进行频谱分析，并画出语音信号的时域与频谱图，并与MATLAB 数字信号处理工具箱中的fft函数进行对比研究，验证自编算法的正确性。
1.2 课程设计要求
1、完成语音信号的采集，利用windows自带的录音机或其他软件，录制一段语音，
时间在1s以内；
2、在MATLAB中编写程序，实现输入信号的倒序；
3、编写程序，实现蝶形运算；
4、画出语音信号的频谱图，与MATLAB数字信号处理工具箱中的fft函数进行对比
研究，并对设计结果进行独立思考和分析；
1.3 课程设计目的
1、学会MATLAB 的使用，掌握MATLAB 的程序设计方法。
2、掌握在Windows 环境下语音信号采集的方法。
3、掌握数字信号处理的基本概念、基本理论和基本方法。
4、掌握MATLAB 设计FIR 和IIR 数字滤波器的方法。
5、学会用MATLAB 对信号进行分析和处理。
1.4 课程设计平台
MATLAB7.1软件
MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动
态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。
2 概要设计
2.1 MATLAB及MATLAB在数字信号领域的应用
MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。附加的工具箱（单独提供的专用MATLAB 函数集）扩展了MATLAB 环境，以解决这些应用领域内特定类型的问题。它以矩阵运算为基础，把计算、可视化、程序设计融合在一个简单易用的交互式工作环境中，是一款数据分析和处理功能都非常强大的工程适用软件。它可以将声音文件变换为离散的数据文件，然后利用其强大的矩阵运算能力处理数据，如数据滤波、傅立叶变换、时域和频域分析、声音回放以及各种图的呈现等，它的信号处理与分析工具箱位语音信号分析提供了十分丰富的功能函数，利用这些功能函数可以快捷而又方便的完成语音信号的处理和分析以及信号的可视化。数字信号处理是MATLAB重要应用的领域之一。
2.2 FFT算法的实现原理
对于有限长序列x(n),若要求其N点的傅里叶变换（DFT）需要经过2N次复数乘法运算和N*(N-1)次复数加法运算。随着N的增加，运算量将急剧增加，而在实际问题中，N往往是较大的，如当N=1024时，完成复数乘法和复数加法的次数分别为百万以上，无论是用通用计算机还是用DSP芯片，都需要消耗大量的时间和机器内存，不能满足实时的要求。因此，DFT的这种运算只能进行理论上的计算，不适合对实时处理要求高的场合。因此，研究作为DSP的快速算法的FFT是相当必要的，快速傅里叶变换（FFT）是为提高DFT运算速度而采用的一种算法，快速算法的种类很多，而且目前仍在改进和提高，它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。基于本学期所学的DIT-FFT的运算规律和编程思想以及MATLAB的
学习和使用，本课设要求在MATLAB环境下编写基2 DIT-FFT算法实现对离散信号的快速傅里叶变换，再与MATLAB软件自带的FFT函数实现对离散信号的傅里叶变换进行比较，如果得到的频谱相同，那么我们编写的程序就是正确的。其中离散信号是通过PC自带的录音机录制一段wav语音信号，用MATLAB采样得到离散序列x1。如果有能力可以选做系统人机对话界面。用GUI界面完成人机交互方便使用的。本课程设计主要是对数字信号的分析。
3 基于MATLAB的FFT算法设计原理
3.1 总体设计流程图
在一个相对较安静的环境下，录下1s左右的wav声音信号，然后对声音进行采样，画出其时域波形和频谱图，其流程图如图1所示：
图1 设计流程图
3.2 语音信号的采集
在实际工作中，我们可以利用windows自带的录音机录制语音文件。采集到语音信号之后，需要对语音信号进行分析，如语音信号的时域分析、频谱分析、语谱图分析。在MATLAB中，我们可以通过[y,fs,bits]=wavread('语音信号路径',[N1 N2])语句。用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。[N1 N2]表示读取从N1点到N2点的值（若只有一个N的点则表示读取前N点的采样值）。向量y 则就代表了一个信号（也即一个复杂的“函数表达式”）也就是说可以像处理一个信号表达式一样处理这个声音信号。
3.3 语音信号的时频分析
利用MATLAB中的“wavread”命令来读入（采集）语音信号，将它赋值给某一向量。再对其进行采样，记住采样频率和采样点数。下面介绍Wavread 函数几种调用格式。
1、y=wavread（file）
功能说明：读取file所规定的wav文件，返回采样值放在向量y中。
2、[y,fs,nbits]=wavread(file)
功能说明：采样值放在向量y中，fs表示采样频率（hz），nbits表示采样位数。3、y=wavread（file，N）
功能说明：读取钱N点的采样值放在向量y中。
4、y=wavread（file，[N1,N2]）
功能说明：读取从N1到N2点的采样值放在向量y中。
接下来，对语音信号speech off.wav进行采样。其程序如下：
[y,fs,nbits]=wavered (…speech off.wav?);
功能说明：把语音信号加载入MATLAB 仿真软件平台中
然后，画出语音信号的时域波形，再对语音信号进行频谱分析。MATLAB提供了快速傅里叶变换算法FFT计算DFT的函数fft,其调用格式如下：
Xk=fft(xn,N)
参数xn为被变换的时域序列向量，N是DFT变换区间长度，当N大于xn的长度时，fft函数自动在xn后面补零。，当N小于xn的长度时，fft函数计算xn的前N个元素，忽略其后面的元素。
原始信号的时域波形图如图2所示：
图2 原始信号的时域波形图原始信号的频域特性图如图3所示：
图3 原始信号的频域特性图
3.4 快速傅里叶变换
快速傅里叶变换（FFT）是为提高DFT运算速度而采用的一种算法。
对一个有限长度序列x(n)的N点的DFT为：X（k）=∑x（n）W^knN (k=0,1,……，N-1；n=0,1,……，N-1;W=e^-j2π/N)
当N=4时，X（k）可展开为：
X（0）= x（0）W^0*4+ x（1）W^0*4 +x（2）W^0*4+ x（3）W^0*4
X（1）= x（0）W^0*4+ x（1）W^1*4 +x（2）W^2*4+ x（3）W^3*4
X（2）= x（0）W^0*4+ x（1）W^2*4 +x（2）W^4*4+ x（3）W^6*4
X（3）= x（0）W^0*4+ x（1）W^3*4 +x（2）W^6*4+ x（3）W^9*4
从上式可以看出,要求4点的DFT,需要16次的复数乘法运算,12次复数乘法运算算。由此类推，要求出N点的DFT，需要N^2次复数乘法运算，N*(N-1)次复数加法运算。当N值较大时，要完成的复数乘法运算和复数加法运算得次数都非常多，无论是用通用计算机还是用DSP芯片，都需要消耗大量的时间，不适合于对实时处理要求高的场合。为了能实时处理DFT,要想减少DFT的运算量可以有两个途径：
第一是降N，N的值减小了，运算量就减少了；第二是利用旋转因子的周期性和对称性，可约性。利用这两个途径实现DFT的快速傅里叶变换（FFT），FFT算法基本上可分为时域抽取法和频域抽取法。
W=e^-j2π/N的性质：
（1）周期性
（2）共轭对称性
（3）可约性
本程序是用基2的按时间抽取的FFT算法（DIT-FFT），设序列x(n)的长度为N，且N满足N=2^M,M为正整数。若N不能满足上述关系，可以将序列x(n)补零实现，则x(n)的N点DFT为：
X（k）=∑x（n）W^knN (k=0,1,……，N-1；n=0,1,……，N-1;W=e^-j2π/N)将n 分为奇数与偶数两部分。
按时间抽取基2-FFT算法的基本思路是将N点序列按时间下标的奇偶分为两个N/2点序列，计算这两个N/2点序列的N/2点DFT，计算量可减小约一半；每一个N/2点序列按照同样的划分原则，可以划分为两个N/4点序列，最后，将原序列划分为多个2点
)
(
)
(N
n
k
N
n
N
k
N
kn
N
W
W
W+
+=
=
*
)
(
*
)
(]
[
]
[n
k
N
n
k
N
kn
N
W
W
W-
-=
=
m
kn
m
N
kn
N
mkn
mN
kn
N
W
W
W
W/
/
,=
=
序列，将计算量大大降低。
按时间下标的奇偶将N点x(n)分别抽取组成两个N/2点序列，分别记为x1(n)和x2(n)，将x(n)的DFT转化为x1(n)和x2(n)的DFT的计算。
利用系数nk
N
W的可约性，即
用蝶形运算可表式为：
以此类推，还可以把x1(n)和x2(n)按n值得奇偶分为两个序列，这样就达到了降N 得目的，从而减少了运算量。
FFT对DFT的数学运算量改进：
直接采用DFT进行计算，运算量为N^2次复数乘法和N*(N-1)次复数乘法。
1
2
,
,1,0
,
)
(
)1
2(
)
(
)
2(
2
1-
=
?
?
?
=
+
=N
r
r
x
r
x
r
x
r
x
()()
()()
()()()
()()()
1
N
21
N N
0,2,4...1,3,5...
11
22
21
2
N N
0,10,1
11
22
21
2
1N2N
0,10,1
221
N
nk
n
N N
nk nk
n n
N N
r k
rk
r r
N N
r k
rk
r r
X k x n W
x n W x n W
x r W x r W
x r W x r W
-
=
--
==
--
+
==
--
+
==
=
=+
=++
=+
∑
∑∑
∑∑
∑∑
2
j
2
j2
22
2
e e
rk
rk
rk rk
N
N N
W W
π
π-
-
===
()()()
11
22
12
22
00
12
,01
N N
rk k rk
N N N
r r
k
N
X k x r W W x r W
X k W X k N
==
=+
=+≤≤-
∑∑
（）（k）
当采用M次FFT时，由N=2^M求得M=logN，运算流图有M级蝶形，每一级都由N/2个蝶形运算构成，这样每一级蝶形运算都需要N/2次复数乘法和N次复数加法。M级运算共需要复数乘法次数为C=N/2*M,复数加法次数为C=N*M。
当N值较大时，FFT减少运算量的特点表现的越明显。
3.4.1 FT的运算规律
3.4.1.1 原位运算
1、N=2^M的FFT共M级运算，每级有N/2蝶形原位计算，当数据输入到存储器以后，每一组蝶形运算后，结果仍然存放在这同一组存储器的同一位置，不需要另辟存储空间，直接最后输出。
2、同一级的蝶形运算每个蝶形运算的输入数据对其他级输入没有影响。
3.4.1.2 倒序运算的规律
输入序列先按自然顺序存入存储单元,然后经变址运算来实现倒位序排列,用J表示倒序的十进制数，对N=2^M，M位的二进制数从左到右各位数权值位N/2，N/4，N/8……2，1。因此，最高位加1相当于J+N/2。
1、如果最高位为0，则直接得到下一个倒序值，J+N/2；
2、最高位为1，则最高位为0(J-N/2)，次高位加1(J+N/4)。
3、类推，直到最后一位二进制数字。
例如，N=8时如下图5所示：
表1 码位倒序（N=8）
倒序的流程图如图4所示：
图4 倒序的流程图
3.4.1.3 规律
设序列x（n）倒序后存放在数组A中，如果蝶形运算两个输入相距B=2^(L-1)。
采用原位运算：蝶形运算可表示为：x(K)=x(K)+x(K+B)*W^Np; x(K+B)=x(K)-x(K+B)*W^Np.其中P=J*2^（M-1）,(J=0,1,2,3……,2^(L-1)-1) 旋转因子的确定：
第L级共有2^(L-1)个旋转因子，以(N=2^3=8为例)：
L=1时，W p
N =W J
L^2
J=0；
L=2时，W p
N =W J
L^2
J=0，1；
L=3时，W p
N =W J
L^2
J=0，1 ，2，3；
W p
N =W J
L^2
=W)
(^2*L
M
J
N
P=J*2^(M-L)，用来确定第L级旋转因子
从输入端开始，共进行M级运算，在进行第L级运算时，依次求出2^(L-1)个旋转因子，然后计算每个旋转因子所对应的2^(M-L)个蝶形元素。第L级的蝶形运算中，每个蝶形运算的两个输入相距B=2^L-1,同一旋转因子对应的蝶形运算相隔2^L个，同一旋转因子对应的蝶形运算有2^M-L个。
注：1、控制第L级顺序运算
2、控制不同种的旋转因子
3、控制同种旋转因子所对应的蝶形运算
蝶形运算的流程图如下图5所示：
图5 蝶形运算的流程图
3.4.2基于MATLAB的FFT所编写程序的框图
基于MATLAB的FFT所编写程序的框图如图5所示：
图6 程序框图
4 调试与结果
自编算法与机带算法仿真波形比较，我们知道MATLAB软件自带FFT算法，我们可以通过比较自编算法仿真结果与机带算法仿真结果来检验自编算法的正确性。
自编算法与FFT算法幅值比较图如图6所示：
图 6 幅值比较图
自编算法FFT算法频谱分贝计较图如图7所示：
图7 分贝比较图
由以上两图可以看出，经过蝶形运算得出的频谱图和信号直接FFT得出的频谱图一致。该程序严格按程序框图编写，思路清晰、容易理解，程序的运行过程在命令窗中一目了然。通过与FFT函数运算的结果比对，程序编写正确。
5 设计总结与心得体会
快速傅氏变换英文名是fast fourier transform，快速傅氏变换（FFT）是离散傅氏变换(DFT)的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。快速傅里叶变换有广泛的应用，如数字信号处理、计算大整数乘法、求解偏微分方程等等。有各种快速算法。
对于复数序列，用离散傅里叶变换。快速傅里叶变换可以计算出与直接计算相同的结果，但只需要的计算复杂度。通常，快速算法要求n能被因数分解，但不是所有的快速傅里叶变换都要求n是合数，对于所有的整数n，都存在快速算法。
除了指数的符号相反、并多了一个1/n的因子，离散傅里叶变换的正变换与逆变换具有相同的形式。因此所有的离散傅里叶变换的快速算法同时适用于正逆变换。
在老师的帮助下我顺利的完成了这个课程设计，通过这次数字信号处理课程设计，让我了解了关于MATLAB软件在数字信号处理方面的应用，又一次学习了MATLAB软件的使用和程序的设计，MATLAB的仿真使我更加深入的了解了数字处理的过程，对我对数字信号处理的理解加深了一步。MATLAB拥有强大的数据仿真能力，在生产和研究中起了非常重要的作用。在这过程中我遇到了所多的难题，通过与老师的交流和学习，让我学会了很多在课堂上没有理解的难点。同时也进一加深了对MATLAB的理解和认识。
为期两周的课设很快接近尾声了，基于MATLAB的FFT算法实现设计已按计划如期全部完成，通过这次DSP课程设计，我对课堂上所学到的理论知识的理解加深了许多，自己动脑、动手设计的能力也得到了较大提高。在这次课程设计的过程中，我对MATLAB 语言有了更深的认识。现在仔细想想，这次课程设计使得我对MATLAB 语言的理解与应用能力得到了较大的提升，也让我认识到只要深入学习，提升的空间永远是存在的。在设计的过程中我遇到了一些问题，如：编写源程序中出现了语法错误等。通过查阅书本和以前设计的程序我发现了产生错误的原因并解决了问题完成了设计。经过反思我发现较大一部分错误是因为操作的不熟练造成的，这也让我明白了要保持设计的高效率必须经常练习。另一方面我也发现了动手实践的重要性。动手实践是理论知识得以灵活运用的必要前提，也是今后走上工作岗位之后能够很好的完成设计工作的技
Matlab中的FFT使用说明
FFT是Fast Fourier Transform（快速傅里叶变换）的简称，FFT算法在MATLAB 中实现的函数是Y=fft（x,n）。刚接触频谱分析用到FFT时，几乎都会对MATLAB 的fft函数产生一些疑惑，下面以看一个例子（根据MATLA帮助修改）。 Fs = 2000; % 设置采样频率 T = 1/Fs; % 得到采用时间 L = 1000; % 设置信号点数，长度1 秒 t = (0:L-1)*T; % 计算离散时间， % 两个正弦波叠加 f1 = 80; A1 = 0.5; % 第一个正弦波100Hz，幅度0.5 f2 = 150; A2 = 1.0 ; % 第2个正弦波150Hz，幅度 1.0 A3 = 0.5; % 白噪声幅度； x = A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t); % 产生离散时间信号； y = x + A3*randn(size(t)); % 叠加噪声； % 时域波形图 subplot(2,1,1) plot(Fs*t(1:50),x(1:50)) title('Sinusoids Signal') xlabel('time (milliseconds)') subplot(2,1,2) plot(Fs*t(1:50),y(1:50)) title('Signal Corrupted with Zero-Mean Random Noise') xlabel('time (milliseconds)') NFFT = 2A nextpow2(L); % 设置FFT点数，一般为2 的N次方，如1024，512 等Y = fft(y,NFFT)/L; % 计算频域信号， f = Fs/2*linspace(0,1,NFFT/2+1); %频率离散化，fft后对应的频率是-Fs/2到Fs/2,由NFFT个离散频点表示 % 这里只画出正频率； % Plot single-sided amplitude spectrum. figure; plot(f,2*abs(Y(1:NFFT/2+1))); % fft 后含幅度和相位，一般观察幅度谱，并把负频率加上去， title('Single-Sided Amplitude Spectrum of y(t)') xlabel('Frequency (Hz)')
实验二 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 级运算的旋转因子
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编译器,从而产生具有交叉列表的汇编程序,而设计者可以对其中的汇编语句进行修改,然后对汇编程序进行编译,产生目标文件。
MATLAB中FFT结果的物理意义
FFT结果的物理意义 最近正在做一个音频处理方面的项目，以前没有学过fft，只是知道有这么个东西，最近这一用才发现原来欠缺这么多，最基本的，连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（ps：横坐标第n个点对应的频率值Fn的计算公式。整个横坐标代表了采样频率Fs，被分为N点。故其频率分辨率为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的结果的模值如图所示。
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 项），直到只剩
用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) 07级电信（2）班 刘坤洋 24 实验一 利用MATLAB 实现信号DFT 的计算 一、实验目的： 1、熟悉利用MATLAB 计算信号DFT 的方法 2、掌握利用MATLAB 实现由DFT 计算线性卷积的方法 二、实验设备：电脑、matlab 软件 三、实验内容： 1、练习用matlab 中提供的内部函数用于计算DFT （1） fft （x ），fft （x ，N ），ifft （x ），ifft （x ，N ）的含义及用法 （2） 在进行DFT 时选取合适的时域样本点数N 请举例，并编程实现 题目： 源程序： >> N=30; %数据的长度 >>L=512; %DFT 的点数 >>f1=100; f2=120; >>fs=600; %抽样频率 >>T=1/fs; %抽样间隔 >>ws=2*pi*fs; >>t=(0:N-1)*T; >>f=cos(4*pi*f1*t)+cos(4*pi*f2*t); >>F=fftshift(fft(f,L)); >>w=(-ws/2+(0:L-1)*ws/L)/(2*pi); >>hd=plot(w,abs(F)); >>ylabel('幅度谱') >> xlabel('频率/Hz') 的频谱 分析利用)π4cos()π4cos()(DFT 21t f t f t x +=Hz 600,Hz 120,Hz 10021===s f f f >> title('my picture') 结果图： （3） 在对信号进行DFT 时选择hamming 窗增加频率分辨率 请举例，并编程实现 题目： 源程序：>> N=50; %数据的长度 >>L=512; %DFT 的点数 >>f1=100;f2=150; >>fs=600; %抽样频率 >>T=1/fs; %抽样间隔 >>ws=2*pi*fs; >>t=(0:N-1)*T; >>f=cos(4*pi*f1*t)+0.15*cos(4*pi*f2*t); 的频谱 分析利用)π4cos(15.0)π4cos()(DFT 21t f t f t x +=Hz 600,Hz 150,Hz 10021===s f f f 实验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.360docs.net/doc/422955922.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.360docs.net/doc/422955922.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 ------＝＋＋＋＋＋ 汇编下的FFT算法实现 ;----------------------------------------------------------- ; 快速付里叶变换子程序 ; ; 入口 : 一维四字节浮点数数组的首地址 ; ; ; 以 2 为基β的值,信号抽样的点数 = 2^β; ; ; 出口 : 在原数组的位置保存结果的值; ; ; 补充说明 : ; ; 1. 该子程序所需 RAM 的容量为 2^β*12 字节,另外需要 ; ; 少量的堆栈作为临时变量的存储空间; ; ; 2. 所需 RAM 空间以输入的首地址为基址,向增加的方向 ; ; 扩展; ; ; 应用举例 : ; ; PUSH #0A000H ; 数组首地址压栈 ; PUSH #6 ; β值压栈 ; CALL FFT ; ;----------------------------------------------------------- PROC FFT FFT: LD FFT_CX,2[SP] LD FFT_AX,#0001H ; SHL FFT_AX,FFT_CL ; 计算采样点数 PUSH FFT_AX ; 将采样点数压栈 SHL FFT_AX,#1 ; LD FFT_BX,6[SP] ; 根据采样点数, ADD FFT_BX,FFT_AX ; 计算出其余 PUSH FFT_BX ; 三个被占用 ADD FFT_BX,FFT_AX ; 空间的首地址; PUSH FFT_BX ; 并依次将首地址 SHL FFT_AX,#1 ; ADD FFT_BX,FFT_AX ; 压栈; PUSH FFT_BX ; LD FFT_CX,6[SP] LD FFT_FX,#2 LOOP1: CLR FFT_AX SUB FFT_EX,FFT_CX,#1 LD FFT_BX,FFT_EX LD FFT_DX,10[SP] LOOP2: 用MATLAB 进行FFT 频谱分析 假设一信号： ()()292.7/2cos 1.0996.2/2sin 1.06.0+++=t t R ππ 画出其频谱图。 分析： 首先，连续周期信号截断对频谱的影响。 DFT 变换频谱泄漏的根本原因是信号的截断。即时域加窗，对应为频域卷积，因此，窗函数的主瓣宽度等就会影响到频谱。 实验表明，连续周期信号截断时持续时间与信号周期呈整数倍关系时，利用DFT 变换可以得到精确的模拟信号频谱。举一个简单的例子： ()ππ2.0100cos +=t Y 其周期为0.02。截断时不同的持续时间影响如图一.1：（对应程序shiyan1ex1.m ） 图 错误!文档中没有指定样式的文字。.1 140.0160.0180.02 截断时，时间间期为周期整数倍，频谱图 0.0250.03 20 40 60 80 100 截断时，时间间期不为周期整数倍，频谱图 其次，采样频率的确定。 根据Shannon 采样定理，采样带限信号采样频率为截止频率的两倍以上，给定信号的采样频率应>1/7.92，取16。 再次，DFT 算法包括时域采样和频域采样两步，频域采样长度M 和时域采样长度N 的关系要符合M ≧N 时，从频谱X(k)才可完全重建原信号。 实验中信号R 经采样后的离散信号不是周期信号，但是它又是一个无限长的信号，因此处理时时域窗函数尽量取得宽一些已接近实际信号。 实验结果如图一.2：其中，0点位置的冲激项为直流分量0.6造成（对应程序为shiyan1.m ） 图 错误!文档中没有指定样式的文字。.2 ?ARMA （Auto Recursive Moving Average ）模型： 将平稳随机信号x(n)看作是零均值，方差为σu 2的白噪声u(n)经过线性非移变系统H(z)后的输出，模型的传递函数为 020406080100120140160180200 0.4 0.50.60.7 0.800.050.10.150.20.250.30.350.40.450.5 50100 150 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 （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); 数字通信课程设计报告书 课题名称 基于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('均方根谱'); 基2的DIT蝶形算法源代码及注释如下： /************FFT***********/ //整个程序输入和输出利用同一个空间x[N]，节约空间 #include 利用MATLAB实现信号DFT的计算
实验2FFT算法实现
(完整word版)stm32F103进行FFT算法教程
128点FFT算法设计方法
汇编语言实现的FFT算法
用MATLAB进行FFT频谱分析
FFT算法实现实验报告
基于matlab的FFT算法程序设计
fft算法代码注释及流程框图