基于matlab的fft算法设计

目录

1引言 (1)

2课程设计内容及要求 (2)

2.1课程设计内容 (2)

2.2课程设计要求 (2)

2.3课程设计目的 (2)

2.3课程设计平台 (2)

3基于MATLAB的FFT算法设计原理 (3)

3.1总体设计流程图 (3)

3.2语音信号的采集 (3)

3.3语音信号的时频分析 (3)

3.4快速傅里叶变换 (6)

3.4.1FFT的运算规律 (8)

3.4.2基于MATLAB的FFT所编写程序的框图 (12)

3.5自编算法与机带算法仿真波形比较 (13)

4设计总结 (16)

参考文献 (17)

附录 (18)

1 引言

随着信息时代，数字时代的到来，数字信号处理已经成为一门极其重要的学科和技术领域。以DSP为核心芯片的处理系统日益变成了数字信号处理系统的主流。它广泛用于电子信息、通信、图像处理、语音处理、生物医学、自动控制、地质探测等领域，受到工程设计和使用人员的青睐。

MATLAB，它是美国Math Works公司推出的一种面向工程和科学计算的交互式计算软件。它以矩阵运算为基础，把计算、可视化、程序设计融合在一个简单易用的交互式工作环境中，是一款数据分析和处理功能都非常强大的工程适用软件。通过本次实习我们学会了分析和处理音频信号，首先要对声音信号进行采集，MATLAB的数据采集工具箱提供了一整套命令和函数，通过调用这些函数和命令，可直接控制声卡进行数据采集。Window自带的录音机程序也可驱动声卡来采集语音信号，并能保存为WA V格式文件，供MATLAB相关函数直接读取、写入或播放。

MATLAB语言是一种数据分析和处理功能十分强大的计算机应用软件，它可以将声音文件变换位离散的数据文件，然后利用其强大的矩阵运算能力处理数据，如数据滤波、傅立叶变换、时域和频域分析、声音回放以及各种图的呈现等，它的信号处理与分析工具箱位语音信号分析提供了十分丰富的功能函数，利用这些功能函数可以快捷而又方便的完成语音信号的处理和分析以及信号的可视化，是人机交互更加便捷。信号处理是MATLAB重要应用的领域之一。

语音信号处理是研究用数字信号处理技术和语音学知识对语音信号进行处理的新兴的学科，是目前发展最为迅速的信息科学研究领域的核心技术之一。通过语音传递信息是人类最重要、最有效、最常用和最方便的交换信息形式。

语音信号的处理与滤波的设计主要是用MATLAB作为工具平台，设计中涉及到声音的录制、播放、存储和读取，语音信号的抽样、频谱分析，滤波器的设计及语音信号的滤波，通过数字信号处理课程的理论知识的综合运用。从实践上初步实现对数字信号的处理。

2 课程设计内容及要求

2.1课程设计内容

录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；在Matlab环境下编写基2 DIT-FFT算法；利用自己编写的算法对已采集的语音信号进行频谱分析，并画出语音信号的时域与频谱图，并与Matlab数字信号处理工具箱中的fft函数进行对比研究，验证自编算法的正确性。

2.2课程设计要求

1.完成语音信号的采集，利用windows自带的录音机或其他软件，录制一段语音，时间在1s以内；

2.在Matlab中编写程序，实现输入信号的倒序；

3.编写程序，实现蝶形运算；

4.画出语音信号的频谱图，与Matlab数字信号处理工具箱中的fft函数进行对

比研究，并对设计结果进行独立思考和分析；

2.3课程设计目的

1.学会MATLAB 的使用，掌握MATLAB 的程序设计方法。

2.掌握在Windows 环境下语音信号采集的方法。

3.掌握数字信号处理的基本概念、基本理论和基本方法。

4.掌握MATLAB 设计FIR 和IIR 数字滤波器的方法。

5.学会用MATLAB 对信号进行分析和处理。

2.4课程设计平台

MATLAB7.1软件

MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。

3 基于MATLAB的FFT算法设计原理

3.1 总体设计流程图

在一个相对较安静的环境下，录下1s左右的wav声音信号，然后对声音进行采样，画出其时域波形和频谱图，其流程图如图1所示：

开始

输入声音信号

对声音信号采样

蝶形运算

原信号fft运算

结束

图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 个元素，忽略其后面的元素。

原始信号的时域波形图如图3所示：

图 3 原始信号的时域波形图原始信号的频域特性图如图4所示：

图 4 原始信号的频域特性图

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点序列，将计算量大大降低。

1. 按时间下标的奇偶将N 点x(n)分别抽取组成两个N/2点序列，分别记为x1(n)和x2(n)，将x(n)的DFT 转化为x1(n)和x2(n)的DFT 的计算。

)

()(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 //,

==12

,

,1,0,

)()12()()2(21-=?

??

=+=N r r x r x r x r x ()()()()()()()()()()1N 0

2

1

N

N

0,2,4...1,3,5 (112)

2

212N

N

0,10,1112

2

2121

N

2N

0,1

0,1

221N 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 -=--==--+==--+====

+

=

++=

+∑∑

∑

∑∑∑∑2j

2j 222

2

e

e

rk N rk

rk rk

N

N N

W W ππ

--===

用蝶形运算可表式为：

以此类推，还可以把x1(n)和x2(n)按n 值得奇偶分为两个序列，这样就达到了降N 得目的，从而减少了运算量。

FFT 对DFT 的数学运算量改进：

直接采用DFT 进行计算，运算量为N^2次复数乘法和N*(N-1)次复数乘法。 当采用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 FFT 的运算规律

（1）原位运算

①N=2^M 的FFT 共M 级运算，每级有N/2蝶形原位计算，当数据输入到存储器以

后，每一组蝶形运算后，结果仍然存放在这同一组存储器的同一位置，不需要另辟存储空间，直接最后输出。

②同一级的蝶形运算每个蝶形运算的输入数据对其他级输入没有影响。

（2）倒序运算的规律

输入序列先按自然顺序存入存储单元,然后经变址运算来实现倒位序排列,用J 表示倒序的十进制数，对N=2^M ，M 位的二进制数从左到右各位数权值位N/2，N/4，N/8……2，1。因此，最高位加1相当于J+N/2。

①.如果最高位为0，则直接得到下一个倒序值，J+N/2； ②.如果最高位为1，则最高位为0(J-N/2)，次高位加1(J+N/4)。 ③.以此类推，直到最后一位二进制数字。

例如 ，N=8时如下图5所示：

()()()1122

122

2

12

,01N

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）

顺序倒序十进制二进制二进制十进制

0 000 000 0

1 001 100 4

2 010 010 2

3 011 110 6

4 100 001 1

5 101 101 5

∑6 110 011 3

7 111 111 7

图5 码位倒序（N=8）

倒序的流程图如图3所示：

图3 倒序的流程图

（3）蝶形运算的运算规律

设序列x （n ）倒序后存放在数组A 中，如果蝶形运算两个输入相距B=2^(L-1).

LH=N/2 N1=N-1 J=0

T=x(J+1) x(J+1)=x(I+1) x(I+1)=T

I=0：N1

I>=J

K=LH

J

J=J+K

J=J-K K=K/2

Y

N

Y

N

采用原位运算：蝶形运算可表示为：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：控制同种旋转因子所对应的蝶形运算

蝶形运算的流程图如图4所示：

：

开始

输入x(n)

N=2^M

倒序

1

L=1:M

B=2^L/2

2

J=0：B-1

P=J*2^(M-L)

3

K=J+1:2^L:N

T=x(K)+x(K+B)* W^Np

x(K+B)=x(K)-x(K+B)* W^Np

x(K)=T

输出

结束

图4 蝶形运算的流程图

3.4.2基于MATLAB的FFT所编写程序的框图

基于MATLAB的FFT所编写程序的框图如图5所示：

图 5 程序框图

3.5 自编算法与机带算法仿真波形比较

我们知道MATLAB 软件自带FFT 算法，我们可以通过比较自编算法仿真结果与

机带算法仿真结果来检验自编算法的正确性。

自编算法与FFT 算法幅值比较图如图6所示：

语音信号采集

完成语音信号时域图

完成语音信号频率特性图

编写fft 程序，画出语音信号频谱图 实现输入信号的倒序 实现一级中不同种蝶形算运 实现一级中相同种蝶形运算

与

m atl ab

自带的

fft

比较

图 6 幅值比较图

自编算法FFT算法频谱分贝计较图如图7所示：

图 7 分贝比较图

由以上两图可以看出，经过蝶形运算得出的频谱图和信号直接FFT得出的频谱图一致。该程序严格按程序框图编写，思路清晰、容易理解，程序的运行过程在命令窗中一目了然。通过与FFT函数运算的结果比对，程序编写正确。

4 设计总结

为期四天的课设很快接近尾声了，基于MATLAB的FFT算法实现设计已按计划如期全部完成，通过这次DSP课程设计，我对课堂上所学到的理论知识的理解加深了许多，自己动脑、动手设计的能力也得到了较大提高。在这次课程设计的过程中，我对 MATLAB 语言有了更深的认识。现在仔细想想，这次课程设计使得我对MATLAB 语言的理解与应用能力得到了较大的提升，也让我认识到只要深入学习，提升的空间永远是存在的。在设计的过程中我遇到了一些问题，如：编写源程序中出现了语法错误等。通过查阅书本和以前设计的程序我发现了产生错误的原因并解决了问题完成了设计。经过反思我发现较大一部分错误是因为操作的不熟练造成的，这也让我明白了要保持设计的高效率必须经常练习。另一方面我也发现了动手实践的重要性。动手实践是理论知识得以灵活运用的必要前提，也是今后走上工作岗位之后能够很好的完成设计工作的技术保证。只有遇到实际问题并根据自己对课堂上获得的专业知识的理解来解决才能真正的提高自己的能力。

这次程序设计让我获益匪浅，对MATLAB也有了进一步的认识：MATLAB功能强、使用灵活方便等。MATLAB是在国内外广泛使用的一种数字信号处理系统，相信除了以上优点，还有许多我还未发现，希望能在以后的学习中有更深入的认识。

在此，我也要感谢一下董翠英老师对我这次课设的大力支持和指导，如果没有老师的帮助我很难在这么短的时间内完成这次课设。

参考文献

[1]范寿康 DSP技术与DSP芯片.北京：电子工业出版社

[2]程佩青.数字信号处理教程.北京：清华大学出版社出版，2001

[3]高西全丁玉美等数字信号处理. 北京：电子工业出版社，2009

[4] 余成波，陶红艳。数字信号处理及MATLAB 实现，北京：清华大学出版社，2008

[5] 曹弋，赵阳。MATLAB 实用教程，北京：电子工业出版社，2007

附录

源程序：

Clear all;

[x,fs,bits]=wavread(‘gl.wav’,2048); %读取声音

x1=reshape(x,1,4096);

sound(x1,fs,bits); %播放语音信号

y1=fft(x1)

N=length(x1); %计数读取信号的点数

t=(1:N)/fs; %信号的时域采样点

figure(1)

plot(t, x1); %画出声音采样后的时域波形

title('原声音信号的时域波形'); %给图形加注标签说明

xlabel('时间/t');

ylabel('振幅/A');

grid ; %添加网格

M=nextpow2(x1); % 求x的长度对应的2的最低幂次m

N=2^M;

if length(x1)

x1=[x1,zeros(1,N-length(x1))]; % 若x的长度不是2的幂，补零到

2的整数幂

end

NV2=N/2;

NM1=N-1;

I=0;

J=0;

while I

if I

T=x1(J+1);

x1(J+1)=x1(I+1);

x1(I+1)=T;

end

K=NV2;

while K<=J

J=J-K;

K=K/2;

end

J=J+K;

I=I+1;

end

%x1;

y=x1; % 将x倒序排列作为y的初始值

WN=exp(-i*2*pi/N); %蝶形运算

for L=1:M

B=2^L/2; %第L级中，每个蝶形的两个输入数据相距

B个点，每级有B个不同的旋转因子for J=0:B-1 % J代表了不同的旋转因子

p=J*2^(M-L);

WNp=WN^p;

for k=J+1:2^L:N % 本次蝶形运算的跨越间隔为2^L

kp=k+B; % 蝶形运算的两个因子对应单元下标的关

系

t=y(kp)*WNp; % 蝶形运算的乘积项

y(kp)=y(k)-t; % 蝶形运算，注意必须先进行减法运算，

然后进行加法运算，否则要使用中间变

量来传递y(k)

y(k)=y(k)+t; % 蝶形运算

end

end

end

%y

figure(2)

[x2,w1]=freqz(x1,1) ; %绘制原始语音信号的频率图

plot(w1/pi,20*log10(abs(x2)));

title('频率特性图')

xlabel('归一化频率');

ylabel('幅度/DB');

grid;

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)')

Simulink下的频谱分析方法及matlab的FFT编程

Simulink下的频谱分析方法 实现功能： 信号发生器一个信号输入，实时显示其频谱分析 调用模块： 信号源（Signal Processing Blockset -> Signal Processing Sources -> Sine Wave） Tip 1：不能用连续的信号源 频谱观察窗（Signal Processing Blockset -> Signal Processing Sources -> Spectrum Scope）Tip 2: 不能用普通的观察窗 Tip 3：必须构上设置中的Buffer input. Buffer size 越大越精细。 Tip 4: 剩下的tips读帮助。 连接关系： 如下图所示 原理框图实验结果：

输出示意图－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 实现功能： 从Workspace读取一组数，进行频谱分析 调用模块： From Workspace Tip 1: 采样时间不能用0，即必须使用离散模式 Tip 2: 从其他模型中Scope保存出来的“Structure with time”的数据可以直接用频谱观察窗（同上一功能） －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 实现功能： 从dSPACE读取一组数，进行频谱分析 实现方法：

1. 从dSPACE读数保存成文件，数据导入Workspace（过程略） 2. 采用从其他模型的Scope保存数据为“Structure with time”的方式构建一个结构变量ScopeData1 3. 使用以下代码将dSPACE数据dscapture拷贝到结构变量ScopeData1中 %% =[0::]; %纯粹为占位，19157为dSPACE保存数据长度 for i=1:19157 end %% 4. 采用下图中的模型进行频谱分析 实验结果： 通过以上方法对单轴压电加速度传感器进行灵敏度分析，下图分别为采用dSPACE和直接利用示波器分析的结果对比。

按时间抽取的基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 级运算的旋转因子

Matlab编程实现FFT变换.

Matlab编程实现FFT变换及频谱分析的程序代码 内容 1．用Matlab产生正弦波,矩形波,以及白噪声信号，并显示各自时域波形图 2．进行FFT变换，显示各自频谱图，其中采样率，频率、数据长度自选 3．做出上述三种信号的均方根图谱,功率图谱,以及对数均方根图谱 4．用IFFT傅立叶反变换恢复信号，并显示恢复的正弦信号时域波形图 源程序 %*************************************************************** **********% % FFT实践及频谱分析% %*************************************************************** **********% %*************************************************************** **********% %***************1.正弦波****************% fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率 %生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(231); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 f=(0:length(y)-1)'*fs/length(y);%进行对应的频率转换 figure(1); subplot(232); plot(f,mag);%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值'); title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱

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的结果的模值如图所示。

MATLAB中FFT的使用方法

MATLAB中FFT的使用方法 一.调用方法 X=FFT(x)； X=FFT(x，N)； x=IFFT(X); x=IFFT(X,N) 用MATLAB进行谱分析时注意： （1）函数FFT返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) →Xk = 39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。 （2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数

n=0:N-1;t=n/fs; %时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %频率序列 subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果：

MATLAB关于FFT频谱分析的程序

MATLAB关于FFT频谱分析的程序 %***************1.正弦波****************% fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率 %生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(231); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 f=(0:length(y)-1)'*fs/length(y);%进行对应的频率转换 figure(1); subplot(232); plot(f,mag);%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值');

title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱 sq=abs(y); figure(1); subplot(233); plot(f,sq); xlabel('频率(Hz)'); ylabel('均方根谱'); title('正弦信号y=2*pi*10t均方根谱'); grid; %求功率谱 power=sq.^2; figure(1); subplot(234); plot(f,power); xlabel('频率(Hz)'); ylabel('功率谱'); title('正弦信号y=2*pi*10t功率谱'); grid; %求对数谱 ln=log(sq); figure(1); subplot(235); plot(f,ln);

利用MATLAB实现信号DFT的计算

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

关于使用Matlab里Powergui的FFTTool分析的问题及解决办法

首先设置 POWERLIB—》powergui，将该模块拖入模型中即可 在需要进行频谱分析的地方连接一示波器 示波器参数设定： Parameters—》Data history—》Save data to workspace； Format—》Structure with time. 运行一次后，双击powergui—》FFT Analysis. 1. 问题1及解决办法 仿真完成后，采用Powergui分析FFT，有时会发生错误："simulation time of the signals is not enough long for the given fundamental frequency". 很多论坛说是仿真时间短了，可能这也是原因，不过更有可能是这样： FFT的数据来自于示波器SCOPE，在SCOPE PARAMETERS/GENERAL选项卡/SAMPLING 中,有DECIMATION和SAMPLE TIME两项，DECIMATION的意思是 The Decimation parameter allows you to write data at every nth sample, where n is the decimation factor. The default decimation, 1, writes data at every time step. 所以，如果选择DECIMATION，记录数据的时刻为第N个采样点，采样点间的时间间隔为采样步长，而在MATLAB Simulink中，如果采用变步长仿真，采样周期就是变化的，这样就很难对采样的数据进行FFT分析，或许软件只认可采样周期一定的数据，所以会出现文首的错误。 如果选择sample time，那么采样周期固定（与仿真步长无关），这样就可以进行FFT 分析了。所以如果遇到文首的错误，可以尝试将示波器的SAMPLing改为sample time，并设定采样周期，Sampling time

实验二 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时：

MATLAB中FFT使用详解

MATLAB中FFT使用详解 一.调用方法 X=FFT(x)； X=FFT(x，N)； x=IFFT(X); x=IFFT(X,N) 用MA TLAB进行谱分析时注意： （1）函数FFT返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) → Xk = 39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。 （2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数 n=0:N-1;t=n/fs; %时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %频率序列 subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅

xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果： fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分：15Hz和40Hz。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析，只需考察0~Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率0~1进行。另外，振幅的大小与所用采样点数有关，采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz

用MATLAB进行FFT频谱分析

用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

MATLAB中FFT的使用方法

MATLAB中FFT的使用方法 调用方法 X=FFT(x)； X=FFT(x,N)；%N为FFT后的数据点数，如果实际信号的数据点数小于N的话，则需要在FFT变换时增加采样点数，或者通过采用频率细分法在原数据后面补充一定数量的0，从而满足N个数据点 X=IFFT(X); X=IFFT(X,N) 一、用MATLAB进行谱分析时注意： （1）函数FFT返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) Xk = 39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。 （2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。

在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。 二、FFT应用举例 例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数 n=0:N-1;t=n/fs; %时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %频率序列 subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100; N=1024;

基于matlab的FFT算法程序设计

数字通信课程设计报告书 课题名称 基于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('均方根谱');

利用MATLAB编写FFT快速傅里叶变换

一、实验目的 1.利用MATLAB 编写FFT 快速傅里叶变换。 2.比较编写的myfft 程序运算结果与MATLAB 中的FFT 的有无误差。 二、实验条件 PC 机，MATLAB7.0 三、实验原理 1. FFT （快速傅里叶变换）原理： 将一个N 点的计算分解为两个N/2点的计算，每个N/2点的计算再进一步分解为N/4点的计算，以此类推。根据DFT 的定义式，将信号x[n]根据采样号n 分解为偶采样点和奇采样点。设偶采样序列为y[n]=x[2n]，奇采样序列为z[n]=x[2n+1]。 上式中的k N W -为旋转因子N k j e /2π-。下式则为y[n]与z[n]的表达式： 2. 蝶形变换的原理：

下图给出了蝶形变换的运算流图，可由两个N/2点的FFT （Y[k]和Z[k]得出N 点FFT X[k]）。同理，每个N/2点的FFT 可以由两个N/4点的FFT 求得。按这种方法，该过程可延迟后推到2点的FFT 。 下图为N=8的分解过程。图中最右边的为8个时域采样点的8点FFTX[k]，由偶编号采样点的4点FFT 和奇编号采样点的4点得到。这4点偶编号又由偶编号的偶采样点的2点FFT 和奇编号的偶采样点的2点FFT 产生。相同的4点奇编号也是如此。依次往左都可以用相同的方法算出，最后由偶编号的奇采样点和奇编号的偶采样点的2点FFT 算出。图中没2点FFT 成为蝶形，第一级需要每组一个蝶形的4组，第二级有每组两个蝶形的两组，最后一级需要一组4个蝶形。 四、实验内容 1.定义函数disbutterfly ，程序根据FFT 的定义：]2 [][][N n x n x n y + +=、n N W N n x n x n z -+ -=])2 [][(][，将序列x 分解为偶采样点y 和奇采样点z 。

MATLAB中FFT的使用方法

MATLAB FFT 的使用方法 2009-08-22 11:00 说明：以下资源来源于《数字信号处理的 MATLAB ；现》万永革主编 一.调用方法 X=FFT(X)； X=FFTX, N); x=IFFT(X); x=IFFT(X,N) 用MATLAB!行谱分析时注意： (1) 函数FFT 返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) Xk = -10.7782 + 6.2929i 7.7071i 5.0000 0 + 5.0000i -10.7782 - 6.2929i Xk 与xn 的维数相同，共有8个元素。Xk 的第一个数对应于直流分量，即频率值 为00 (2) 做FFT 分析时，幅值大小与FFT 选择的点数有关，但不影响分析结果。在 IFFT 时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果 乘以2除以N 即可。 二.FFT 应用举例 例 1: x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t) 。采样频率 fs=100Hz,分别绘 制N=128 1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数 n=0:N-1;t=n/fs; %时间序歹 U x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); % 信号 39.0000 5.0000i 4.7782 + 7.7071i 4.7782

y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %? 率序列 subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N⑵,mag(1:N⑵)；％绘出Nyquist 频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; 以对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); % 信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag); % 绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2)); % 绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果:

用matlab进行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的结果的模值如图所示。