离散傅里叶变换和快速傅里叶变换

实验报告课程名称： 信号分析与处理 指导老师： 成绩：__________________实验名称：离散傅里叶变换和快速傅里叶变换 实验类型： 基础实验 同组学生姓名：第二次实验 离散傅里叶变换和快速傅里叶变换一、实验目的1.1掌握离散傅里叶变换（DFT ）的原理和实现；1.2掌握快速傅里叶变换（FFT ）的原理和实现，掌握用FFT 对连续信号和离散信号进行谱分析的方法。

1.3 会用Matlab 软件进行以上练习。

二、实验原理2.1关于DFT 的相关知识序列x (n )的离散事件傅里叶变换（DTFT ）表示为nj n j en x e X Ω-∞-∞=Ω∑=)()(，如果x (n )为因果有限长序列，n =0,1,...,N-1，则x (n )的DTFT 表示为n j N n j e n x e X Ω--=Ω∑=1)()(，x (n )的离散傅里叶变换（DFT ）表达式为)1,...,1,0()()(21-==--=∑N k en x k X nk NjN n π，序列的N 点DFT 是序列DTFT 在频率区间[0,2π]上的N 点灯间隔采样，采样间隔为2π/N 。

通过DFT ，可以完成由一组有限个信号采样值x (n )直接计算得到一组有限个频谱采样值X (k )。

X (k )的幅度谱为)()()(22k X k X k X I R +=，其中下标R 和I 分别表示取实部、虚部的运算。

X (k )的相位谱为)()(arctan)(k X k X k R I =ϕ。

离散傅里叶反变换（IDFT ）定义为)1,...,1,0()(1)(21-==∑-=N n e k X N n x nk Nj N n π。

2.2关于FFT 的相关知识快速傅里叶变换（FFT ）是DFT 的快速算法，并不是一个新的映射。

FFT 利用了n Nj eπ2-函数的周期性和对称性以及一些特殊值来减少DFT 的运算量，可使DFT 的运算量下降几个数量级，从而使数字信号处装订线理的速度大大提高。

若信号是连续信号，用FFT 进行谱分析时，首先必须对信号进行采样，使之变成离散信号，然后就可以用FFT 来对连续信号进行谱分析。

为了满足采样定理，一般在采样之前要设置一个抗混叠低通滤波器，且抗混叠滤波器的截止频率不得高于与采样频率的一半。

比较DFT 和IDFT 的定义，两者的区别仅在于指数因子的指数部分的符号差异和幅度尺度变换，因此可用FFT 算法来计算IDFT 。

三、实验内容与相关分析（共6道）说明：为了便于老师查看，现将各题的内容写在这里——题目按照3.1、3.2、...、3.6排列。

每道题包含如下内容：题干、解答（思路、M 文件源代码、命令窗口中的运行及其结果）、分析。

其中“命令窗口中的运行及其结果”按照小题顺序排列，各小题包含命令与结果（图形或者序列）。

3.1 求有限长离散时间信号x (n )的离散时间．．傅里叶变换（DTFT ）X (e j Ω)并绘图。

（1）已知⎩⎨⎧≤≤-=其他0221)(n n x ；（2）已知1002)(≤≤=n n x n。

【解答】思路：这是DTFT 的变换，按照定义编写DTFT 的M 文件即可。

考虑到自变量Ω是连续的，为了方便计算机计算，计算时只取三个周期[-2π,4π]中均匀的1000个点用于绘图。

理论计算的各序列DTFT 表达式，请见本题的分析。

M 文件源代码（我的Matlab 源文件不支持中文注释，抱歉）： function DTFT(n1,n2,x)%This is a DTFT function for my experiment of Signal Processing & Analysis. w=0:2*pi/1000:2*pi;%Define the bracket of omega for plotting. X=zeros(size(w));%Define the initial values of X. for i=n1:n2X=X+x(i-n1+1)*exp((-1)*j*w*i);%It is the definition of DTFT. endAmp=abs(X);%Acquire the amplification.Phs=angle(X);%Acquire the phase angle (radian). subplot(1,2,1);plot(w,Amp,'r'); xlabel('\Omega');ylabel('Amplification');hold on ; %Plot amplification on the left. subplot(1,2,2);plot(w,Phs,'b');xlabel('\Omega');ylabel('Phase Angle (radian)');hold off ; %Plot phase angle on the right. end命令窗口中的运行及其结果（理论计算的各序列DTFT 表达式，请见本题的分析）： 第（1）小题>> n=(-2:2); >> x=1.^n;>> DTFT(-2,2,x);-5051000.511.522.533.544.55ΩA m p l i f i c a t i on-50510-4-3-2-101234ΩP h a s e A n g l e (r a d i a n )第（2）小题>> n=(0:10); >> x=2.^n;>> DTFT(0,10,x);-55106008001000120014001600180020002200ΩA m p l i f i c a t i on-50510-4-3-2-101234ΩP h a s e A n g l e (r a d i a n )【分析】对于第（1）小题，由于序列x(n)只在有限区间（-2，-1，-，1，2）上为1，所以是离散非周期的信号。

它的幅度频谱相应地应该是周期连续信号。

事实上，我们可计算出它的表达式：()Ω-Ω-Ω-Ω-Ω-=Ω-∞+-∞=Ω---=Ω⇒--===Ω∑∑j j j j j n nj n n j X n x X e1e 1)(e 1e 1e e e )()(55222，可见幅度频谱拥有主极大图3.1.1在[-2π,4π]范围内3个周期的幅度谱和相位谱（弧度制）图3.1.2在[-2π,4π]范围内3个周期的幅度谱和相位谱（弧度制）和次极大，两个主极大间有|5-1|=4个极小，即有3个次级大。

而对于它的相位频谱，则是周期性地在-π、0、π之间震荡。

对于第（2）小题，由于是离散非周期的信号。

它的幅度频谱相应地应该是周期连续信号。

而它的表达式：()Ω-Ω-Ω-=Ω-+∞-∞=Ω--≈Ω⇒--===Ω∑∑j j j n nj n nj X n x X e212)(e 21e 212ee)()(11111110，因此主极大之间只有|0-1|=1个极小，不存在次级大。

而对于它的相位频谱，则是在一个长为2π的周期内有|11-1|=10次振荡。

而由DTFT 的定义可知，频谱都是以2π为周期向两边无限延伸的。

由于DTFT 是连续谱，对于计算机处理来说特别困难，因此我们才需要离散信号的频谱也离散，由此构造出DFT （以及为加速计算DFT 的FFT ）。

3.2已知有限长序列x (n )={8,7,9,5,1,7,9,5},试分别采用DFT 和FFT 求其离散傅里叶变换X (k )的幅度、相位图。

【解答】思路：按照定义编写M 文件即可。

M 文件源代码： i) DFT 函数：function DFT(N,x)%This is a DFT function for my experiment of Signal Processing & Analysis. k=(0:N-1);%Define variable k for DFT.X=zeros(size(k));%Define the initial valves of X. for i=0:N-1X=X+x(i+1)*exp((-1)*j*2*k*pi/N*i);%It is the definition of DFT. endAmp=abs(X);%Acquire the amplification.Phs=angle(X);%Acquire the phase angle (radian). subplot(1,2,1);stem(k,Amp,'.',’MarkerSize ’,18); xlabel('k');ylabel('Amplification');hold on ; %Plot amplification on the left. subplot(1,2,2);stem(k,Phs,'*');xlabel('k');ylabel('Phase Angle (radian)');hold off ; %Plot phase angle on the right. endii) 基2-FFT 函数function myFFT(N,x)%This is a base-2 FFT function. lov=(0:N-1); j1=0;for i=1:N %indexed addressing if i<j1+1temp=x(j1+1); x(j1+1)=x(i); x(i)=temp; endk=N/2;while k<=j1j1=j1-k;k=k/2;endj1=j1+k;enddigit=0;k=N;while k>1digit=digit+1;k=k/2;endn=N/2;% Now we start the "butterfly-shaped" process.for mu=1:digitdif=2^(mu-1);%Differnce between the indexes of the target variables.idx=1;for i=1:nidx1=idx;idx2=1;for j1=1:N/(2*n)r=(idx2-1)*2^(digit-mu);wn=exp(j*(-2)*pi*r/N);%It is the "circulating coefficients".temp=x(idx);x(idx)=temp+x(idx+dif)*wn;x(idx+dif)=temp-x(idx+dif)*wn;idx=idx+1;idx2=idx2+1;endidx=idx1+2*dif;endn=n/2;endAmp=abs(x);%Acquire the amplification.Phs=angle(x);%Acquire the phase angle (radian).subplot(1,2,1);stem(lov,Amp,'.',’MarkerSize’,18);xlabel('FFT k');ylabel('FFT Amplification');hold on;%Plot the amplification.subplot(1,2,2);stem(lov,Phs,'*');xlabel('FFT k');ylabel('FFT Phase Angle (radian)');hold off; end命令窗口中的运行及其结果：DFT ：>> x=[8,7,9,5,1,7,9,5]; >> DFT(8,x);24680102030405060kA m p l i f i c a t i o n02468-3-2-1123kP h a s e A n g l e (r a d i a n )FFT ： >> x=[8,7,9,5,1,7,9,5]; >> myFFT(8,x);24680102030405060FFT kF F T A m p l i f i c a t i o n02468-3-2-1123FFT kF F T P h a s e A n g l e (r a d i a n )【分析】DFT 是离散信号、离散频谱之间的映射。