离散傅里叶变换和快速傅里叶变换.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
ifitemp=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
升级会员 