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

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

《离散傅里叶变换和快速傅里叶变换.docx》由会员分享，可在线阅读，更多相关《离散傅里叶变换和快速傅里叶变换.docx（33页珍藏版）》请在冰豆网上搜索。

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

实验报告

课程名称：

信号分析与处理指导老师：

成绩：

__________________

实验名称：

离散傅里叶变换和快速傅里叶变换实验类型：

基础实验同组学生姓名：

第二次实验离散傅里叶变换和快速傅里叶变换

一、实验目的

1.1掌握离散傅里叶变换（DFT）的原理和实现；

1.2掌握快速傅里叶变换（FFT）的原理和实现，掌握用FFT对连续信号和离散信号进行谱分析的方法。

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

二、实验原理

2.1关于DFT的相关知识

序列x（n）的离散事件傅里叶变换（DTFT）表示为

，

如果x（n）为因果有限长序列，n=0,1,...,N-1，则x（n）的DTFT表示为

，

x（n）的离散傅里叶变换（DFT）表达式为

，

序列的N点DFT是序列DTFT在频率区间[0,2π]上的N点灯间隔采样，采样间隔为2π/N。

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

X（k）的幅度谱为

，其中下标R和I分别表示取实部、虚部的运算。

X（k）的相位谱为

。

离散傅里叶反变换（IDFT）定义为

。

2.2关于FFT的相关知识

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

FFT利用了

函数的周期性和对称性以及一些特殊值来减少DFT的运算量，可使DFT的运算量下降几个数量级，从而使数字信号处理的速度大大提高。

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

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

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

3、实验内容与相关分析（共6道）

说明：

为了便于老师查看，现将各题的内容写在这里——

题目按照3.1、3.2、...、3.6排列。

每道题包含如下内容：

题干、解答（思路、M文件源代码、命令窗口中的运行及其结果）、分析。

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

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

（1）已知

；

（2）已知

。

【解答】

思路：

这是DTFT的变换，按照定义编写DTFT的M文件即可。

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

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

M文件源代码（我的Matlab源文件不支持中文注释，抱歉）：

functionDTFT（n1,n2,x）

%ThisisaDTFTfunctionformyexperimentofSignalProcessing&Analysis.

w=0:

2*pi/1000:

2*pi;%Definethebracketofomegaforplotting.

X=zeros（size（w））;%DefinetheinitialvaluesofX.

fori=n1:

n2

X=X+x（i-n1+1）*exp（（-1）*j*w*i）;%ItisthedefinitionofDTFT.

end

Amp=abs（X）;%Acquiretheamplification.

Phs=angle（X）;%Acquirethephaseangle（radian）.

subplot（1,2,1）;

plot（w,Amp,'r'）;xlabel（'\Omega'）;ylabel（'Amplification'）;holdon;

%Plotamplificationontheleft.

subplot（1,2,2）;

plot（w,Phs,'b'）;xlabel（'\Omega'）;ylabel（'PhaseAngle（radian）'）;holdoff;

%Plotphaseangleontheright.

end

命令窗口中的运行及其结果（理论计算的各序列DTFT表达式，请见本题的分析）：

第

（1）小题

>>n=（-2:

2）;

>>x=1.^n;

>>DTFT（-2,2,x）;

图3.1.1在[-2π,4π]范围内3个周期的幅度谱和相位谱（弧度制）

第

（2）小题

>>n=（0:

10）;

>>x=2.^n;

>>DTFT（0,10,x）;

图3.1.2在[-2π,4π]范围内3个周期的幅度谱和相位谱（弧度制）

【分析】

对于第

（1）小题，由于序列x（n）只在有限区间（-2，-1，-，1，2）上为1，所以是离散非周期的信号。

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

事实上，我们可计算出它的表达式：

，可见幅度频谱拥有主极大和次极大，两个主极大间有|5-1|=4个极小，即有3个次级大。

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

对于第

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

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

而它的表达式：

，因此主极大之间只有|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函数：

functionDFT（N,x）

%ThisisaDFTfunctionformyexperimentofSignalProcessing&Analysis.

k=（0:

N-1）;%DefinevariablekforDFT.

X=zeros（size（k））;%DefinetheinitialvalvesofX.

fori=0:

N-1

X=X+x（i+1）*exp（（-1）*j*2*k*pi/N*i）;%ItisthedefinitionofDFT.

end

Amp=abs（X）;%Acquiretheamplification.

Phs=angle（X）;%Acquirethephaseangle（radian）.

subplot（1,2,1）;

stem（k,Amp,'.',’MarkerSize’,18）;xlabel（'k'）;ylabel（'Amplification'）;holdon;

%Plotamplificationontheleft.

subplot（1,2,2）;

stem（k,Phs,'*'）;xlabel（'k'）;ylabel（'PhaseAngle（radian）'）;holdoff;

%Plotphaseangleontheright.

end

ii）基2-FFT函数

functionmyFFT（N,x）

%Thisisabase-2FFTfunction.

lov=（0:

N-1）;

j1=0;

fori=1:

N%indexedaddressing

ifi

temp=x（j1+1）;

x（j1+1）=x（i）;

x（i）=temp;

end

k=N/2;

whilek<=j1

j1=j1-k;

k=k/2;

end

j1=j1+k;

end

digit=0;

k=N;

whilek>1

digit=digit+1;

k=k/2;

end

n=N/2;%Nowwestartthe"butterfly-shaped"process.

formu=1:

digit

dif=2^（mu-1）;%Differncebetweentheindexesofthetargetvariables.

idx=1;

fori=1:

n

idx1=idx;

idx2=1;

forj1=1:

N/（2*n）

r=（idx2-1）*2^（digit-mu）;

wn=exp（j*（-2）*pi*r/N）;%Itisthe"circulatingcoefficients".

temp=x（idx）;

x（idx）=temp+x（idx+dif）*wn;

x（idx+dif）=temp-x（idx+dif）*wn;

idx=idx+1;

idx2=idx2+1;

end

idx=idx1+2*dif;

end

n=n/2;

end

Amp=abs（x）;%Acquiretheamplification.

Phs=angle（x）;%Acquirethephaseangle（radian）.

subplot（1,2,1）;

stem（lov,Amp,'.',’MarkerSize’,18）;

xlabel（'FFTk'）;ylabel（'FFTAmplification'）;holdon;

%Plottheamplification.

subplot（1,2,2）;

stem（lov,Phs,'*'）;xlabel（'FFTk'）;ylabel（'FFTPhaseAngle（radian）'）;holdoff;

end

命令窗口中的运行及其结果：

DFT：

>>x=[8,7,9,5,1,7,9,5];

>>DFT（8,x）;

图3.2.1DFT的幅度谱和相位谱（弧度制）

FFT：

>>x=[8,7,9,5,1,7,9,5];

>>myFFT（8,x）;

图3.2.2FFT算法的幅度谱和相位谱（弧度制）

图3.2.1DFT的幅度谱和相位谱（相位是弧度制的）

【分析】

DFT是离散信号、离散频谱之间的映射。

在这里我们可以看到序列的频谱也被离散化。

事实上，我们可以循着DFT构造的方法验证这个频谱：

首先，将序列做N=8周期延拓，成为离散周期信号。

然后利用DFS计算得到延拓后的频谱：

，从而取DFS的主值区间得到DFT，与图一致。

因此计算正确。

而对于FFT，我们可以看到它给出和DFT一样的结果，说明了FFT算法就是DFT的一个等价形式。

不过，由于序列不够长，FFT在计算速度上的优越性尚未凸显。

3.3已知连续时间信号x（t）=3cos8πt,X（ω）=

，该信号从t=0开始以采样周期Ts=0.1s进行采样得到序列x（n），试选择合适的采样点数，分别采用DFT和FFT求其离散傅里叶变换X（k）的幅度、相位图，并将结果与X（k）的幅度、相位图，并将结果与X（ω）相比较。

【解答】

思路：

此题与下一题都是一样的操作，可以在编程时统一用变量g（0或1）来控制是否有白噪声。

这里取g=0（无白噪声）。

另外，分别取12点、20点、28点采样，以考察采样长度的选择与频谱是否泄漏的关系。

M文件源代码：

i）采样函数：

functionxs=sampJune3（N,Ts,g）

%ThisisafunctionappliedtoProblem3&4.

%N:

numberofsamplingpoints;Ts:

samplingperiod;g=0:

Nogaussiannoise;g=1:

gussiannoiseexists.

n=1:

N;

fori=1:

N%Notethatimuststartat0inanalysis.ThusImadeaadaptation.

sample（i）=3*cos（8*pi*Ts*（i-1））+g*randn;

%InMatlab,indexstartsat1,whichisnotconsistentwithourhabit.ThusImadeaadaptationincodes.

%Itisasamplingprocesswith（g=1）/without（g=0）noise.

end

xs=sample;

plot（n-1,sample,'.',’MarkerSize’,18）;xlabel（'n'）;ylabel（'value'）;holdoff;

%Mustuse（n-1）,becauseinMatlab,indexstartsat1.

end

ii）DFT和基2-FFT函数的代码，请见第3.2节。

不需再新编一个。

命令窗口中的运行及其结果：

12点采样：

>>xs=sampJune3（12,0.1,0）;%末尾的0表示无噪声。

>>DFT（12,xs）;

>>myFFT（12,xs）;

图3.3.1进行12点采样得到的序列

图3.3.2DFT的幅度谱和相位谱（弧度制），出现了泄漏

图3.3.3FFT的幅度谱和相位谱（弧度制）。

出现了频谱泄漏。

20点采样：

>>xs=sampJune3（20,0.1,0）;%末尾的0表示无噪声。

>>DFT（20,xs）;

>>myFFT（20,xs）;

图3.3.4进行20点采样得到的序列

图3.3.5DFT的幅度谱和相位谱（弧度制）。

频谱无泄漏。

图3.3.6FFT的幅度谱和相位谱（弧度制）。

频谱无泄漏。

28点采样：

>>xs=sampJune3（28,0.1,0）;%末尾的0表示无噪声。

>>DFT（28,xs）;

>>myFFT（28,xs）;

图3.3.7进行28点采样得到的序列

图3.3.8DFT的幅度谱和相位谱（弧度制）。

再次出现频谱泄漏。

图3.3.9FFT的幅度谱和相位谱（弧度制）。

再次出现频谱泄漏。

【分析】

分别取12点、20点、28点采样，以考察采样长度的选择与频谱是否泄漏之间的关系。

现在与原信号频谱

比较后可以得出如下结论：

图3.3.10原信号的频谱（由两个冲激函数组成）

原信号的频谱是

，在±8π上各有一强度为3π的谱线，在其余频率上为0。

可见原信号被0.1s采样周期的采样信号离散化之后，谱线以20π为周期重复，并且只在（20k±8）π（k为整数）处非0。

那么，在20点DFT（采样时间原信号周期的整数倍）中，只有第8根、第12根谱线非0。

而在12点、28点DFT中，由于采样时间不是原信号周期的整数倍，谱线将向两边泄漏。

不过，对比12点采样和28点采样，我们还可以发现，28点采样频谱的主谱线高度是次谱线高度的4倍，儿12点采样频谱的主谱线高度是次谱线高度的3倍。

可见，在无法保证采样时间是信号周期整数倍的情况下，增加采样时间有助于减轻频谱泄漏的程度。

3.4对第3步中所述连续时间信号叠加高斯白噪声信号，重复第3步过程。

【解答】

思路：

此题与上一题都是一样的操作，可以在编程时统一用变量g（0或1）来控制是否有白噪声。

这里取g=1（有白噪声）。

另外，仍然分别取12点、20点、28点采样，以考察采样长度的选择与频谱是否泄漏的关系。

M文件源代码：

不需要再新编程序。

可以直接引用上面的函数：

sampJune3（N,Ts,g），取g=1，以体现存在白噪声

DFT（N,x）

myFFT（N,x）

命令窗口中的运行及其结果：

12点采样：

>>xs=sampJune3（12,0.1,1）;%末尾的1表示有噪声。

>>DFT（12,xs）;

>>myFFT（12,xs）;

图3.4.1进行12点采样得到的含噪声的序列

图3.4.2含噪声序列DFT的幅度谱和相位谱（弧度制）。

图3.4.3含噪声FFT的幅度谱和相位谱（弧度制）。

20点采样：

>>xs=sampJune3（20,0.1,1）;%末尾的1表示有噪声。

>>DFT（20,xs）;

>>myFFT（20,xs）;

图3.4.4进行20点采样得到的含噪声序列

图3.4.5含噪声DFT的幅度谱和相位谱（弧度制）。

图3.4.6含噪声FFT的幅度谱和相位谱（弧度制）。

28点采样：

>>xs=sampJune3（28,0.1,0）;%末尾的1表示有噪声。

>>DFT（28,xs）;

>>myFFT（28,xs）;

图3.4.7进行28点采样得到的含噪声序列

图3.4.8含噪声DFT的幅度谱和相位谱（弧度制）。

图3.4.9含噪声FFT的幅度谱和相位谱（弧度制）。

【分析】

依然分别取12点、20点、28点采样。

仍然与原信号的频谱

（图3.3.10）比较，可以得到结论：

由于叠加了噪声，所以频谱都受到了一定的干扰。

由于白噪声在各个频率的功率相等，因此频谱上各处的干扰也是均匀随机的。

不过，通过对比我们可以发现，20点采样（无噪声时不发生泄漏的采样方法）在存在噪声时，仍然可以明显区分出原信号的谱线。

第二好的是28点采样，因为采样时间较长，即使存在频谱泄漏也能较好地区分原信号的谱线。

而最差的是12点采样，由于噪声的存在和严重的频谱泄漏，它的次谱线与主谱线的高度相差不大，使原信号不明显。

3.5已知序列

，X（k）是x（n）的6点DFT，设

。

（1）若有限长序列y（n）的6点DFT是

，求y（n）。

（2）若有限长序列w（n）的6点DFTW（k）是

的实部，求w（n）。

（3）若有限长序列q（n）的3点DFT是

，k=0,1,2，求q（n）。

【解答】

思路：

这是对DFT进行变换后求IDFT。

考虑到IDFT和DFT定义的对称性，可以在DFT的基础上略加调整既可用于计算。

首先，∵

，

∴它的6点采样是序列是

。

值得指出的是，在Matlab中，数组的序号是从1开始的（而在信号分析中习惯从0开始），不过我在上面编程时已考虑到这一情况，具体可见实验报告最后的“附录”。

首先生成x（n）的6点DFT，再按照各小题分别转换，最后求相应的IDFT。

M文件源代码：

i）输出x（n）的6点DFT的函数：

functionX=ExportDFT（N,x）

%ThisisaDFTfunctionthatexportsthesolutiontovectorX.

k=（0:

N-1）;%DefinevariablekforDFT.

X=zeros（size（k））;%DefinetheinitialvalvesofX.

fori=0:

N-1

X=X+x（i+1）*exp（（-1）*j*2*k*pi/N*i）;%ItisthedefinitionofDFT.

end

end

ii）算第

（1）小题的Y（k）的函数：

functionY=Convertor1（X）

%Thisisamathematicalconvertorforthesubproblem

（1）.

fork=1:

6

Y（k）=exp（（-j）*2*pi*（-4*（k-1））/6）*X（k）;

%Herewemustuse（k-1）,becauseinMatlabindexstartsat1.

end

end

iii）算第

（2）小题的W（k）的函数：

functionW=Convertor2（X）

%Thisisamathematicalconvertorforthesubproblem

（2）.

W=real（X）;%AcquiretherealpartofX.

end

iv）算第（3）小题的Q（k）的函数：

functionQ=Convertor3（X）

%Thisisamathematicalconvertorforthesubproblem（3）.

Q=zeros（3）;%Detailedexplanationgoeshere

fortmp=1:

3

Q（tmp）=X（2*tmp）;

end

end

v）输出IDFT的函数：

functionx=ExportIDFT（N,X）

%ThisistheIDFTfunctionformyexperiment.

n=（0:

N-1）;%DefinevariablenforIDFT.

x=zeros（size（n））;%Definetheinitialvalvesofx.

fork=0:

N-1

x=x+X（k+1）*exp（j*2*k*pi/N*n）;

end

x=x/N;

a=real（x）;%WeMUSTusereal（x）,thoughwemayALREADYknowxisreal.

%It'scausedbynumericcalculation（notanalyticcalculation）inMatlab.

stem（n,a,'.','MarkerSize',18）;xlabel（'n'）;ylabel（'Solution'）;

end

命令窗口中的运行及其结果：

>>x=[4,3,2,1,0,0];

>>X=ExportDFT（6,x）;

第

（1）小题

>>Y=Convertor1（X）;

>>y=ExportIDFT（6,Y）

y=

Columns1through4

0.0000+0.0000i0.0000+0.0000i4.0000+0.0000i3.0000+0.0000i

%虚部都是0，说明是实数

Columns5through6

2.0000+0.0000i1.0000-0.0000i%虚部都是0，说明是实数

%事实上，在Matlab中，由于数值计算的截断误差，对原复数做乘法后，答案的虚部可能有一极小的量。

答案：

y（n）={0,0,4,3,2,1}

图3.5.1输出的y（n），这是对x（n）的圆周移位。

第

（2）小题

>>W=Convertor2（X）;

>>w=ExportIDFT（6,W）

w=

Columns1through4

4.0000+0.0000i1.5000+0.0000i1.0000+0.0000i1.0000+0.0000i

%虚部都是0，说明是实数

Columns5through6

1.0000+0.0000i1.5000+0.0000i%虚部都是0，说明是实数；

%事实上，在Matlab中，由于数值计算的截断误差，对原复数做乘法后，答案的虚部可能有一极小的量。

答案：

w（n）={0,0,4,3,2,1}

图3.5.2输出的w（n