离散傅里叶变换和快速傅里叶变换Word下载.docx
《离散傅里叶变换和快速傅里叶变换Word下载.docx》由会员分享，可在线阅读，更多相关《离散傅里叶变换和快速傅里叶变换Word下载.docx（33页珍藏版）》请在冰点文库上搜索。
题目按照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（'
PhaseAngle（radian）'
holdoff;
%Plotphaseangleontheright.
命令窗口中的运行及其结果（理论计算的各序列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&
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.
stem（k,Amp,'
.'
’MarkerSize’,18）;
k'
stem（k,Phs,'
*'
ii）基2-FFT函数
functionmyFFT（N,x）
%Thisisabase-2FFTfunction.
lov=（0:
j1=0;
fori=1:
N%indexedaddressing
ifi<
j1+1
temp=x（j1+1）;
x（j1+1）=x（i）;
x（i）=temp;
end
k=N/2;
whilek<
=j1
j1=j1-k;
k=k/2;
j1=j1+k;
digit=0;
k=N;
whilek>
1
digit=digit+1;
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;
idx=idx1+2*dif;
n=n/2;
Amp=abs（x）;
Phs=angle（x）;
stem（lov,Amp,'
FFTk'
FFTAmplification'
%Plottheamplification.
stem（lov,Phs,'
FFTPhaseAngle（radian）'
命令窗口中的运行及其结果：
DFT：
x=[8,7,9,5,1,7,9,5];
DFT（8,x）;
图3.2.1DFT的幅度谱和相位谱（弧度制）
FFT：
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点采样，以考察采样长度的选择与频谱是否泄漏的关系。
i）采样函数：
functionxs=sampJune3（N,Ts,g）
%ThisisafunctionappliedtoProblem3&
4.
%N:
numberofsamplingpoints;
Ts:
samplingperiod;
g=0:
Nogaussiannoise;
g=1:
gussiannoiseexists.
n=1:
N;
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.
xs=sample;
plot（n-1,sample,'
n'
value'
%Mustuse（n-1）,becauseinMatlab,indexstartsat1.
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）;
DFT（20,xs）;
myFFT（20,xs）;
图3.3.4进行20点采样得到的序列
图3.3.5DFT的幅度谱和相位谱（弧度制）。
频谱无泄漏。
图3.3.6FFT的幅度谱和相位谱（弧度制）。
28点采样：
xs=sampJune3（28,0.1,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点采样，以考察采样长度的选择与频谱是否泄漏的关系。
不需要再新编程序。
可以直接引用上面的函数：
sampJune3（N,Ts,g），取g=1，以体现存在白噪声
DFT（N,x）
myFFT（N,x）
xs=sampJune3（12,0.1,1）;
%末尾的1表示有噪声。
图3.4.1进行12点采样得到的含噪声的序列
图3.4.2含噪声序列DFT的幅度谱和相位谱（弧度制）。
图3.4.3含噪声FFT的幅度谱和相位谱（弧度制）。
xs=sampJune3（20,0.1,1）;
图3.4.4进行20点采样得到的含噪声序列
图3.4.5含噪声DFT的幅度谱和相位谱（弧度制）。
图3.4.6含噪声FFT的幅度谱和相位谱（弧度制）。
图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。
i）输出x（n）的6点DFT的函数：
functionX=ExportDFT（N,x）
%ThisisaDFTfunctionthatexportsthesolutiontovectorX.
%DefinevariablekforDFT.
%DefinetheinitialvalvesofX.
%ItisthedefinitionofDFT.
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.
iii）算第
（2）小题的W（k）的函数：
functionW=Convertor2（X）
%Thisisamathematicalconvertorforthesubproblem
（2）.
W=real（X）;
%AcquiretherealpartofX.
iv）算第（3）小题的Q（k）的函数：
functionQ=Convertor3（X）
%Thisisamathematicalconvertorforthesubproblem（3）.
Q=zeros（3）;
%Detailedexplanationgoeshere
fortmp=1:
3
Q（tmp）=X（2*tmp）;
v）输出IDFT的函数：
functionx=ExportIDFT（N,X）
%ThisistheIDFTfunctionformyexperiment.
n=（0:
%DefinevariablenforIDFT.
x=zeros（size（n））;
%Definetheinitialvalvesofx.
fork=0:
x=x+X（k+1）*exp（j*2*k*pi/N*n）;
x=x/N;
a=real（x）;
%WeMUSTusereal（x）,thoughwemayALREADYknowxisreal.
%It'
scausedbynumericcalculation（notanalyticcalculation）inMatlab.
stem（n,a,'
'
MarkerSize'
18）;
Solution'
x=[4,3,2,1,0,0];
X=ExportDFT（6,x）;
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）的圆周移位。
W=Convertor2（X）;
w=ExportIDFT（6,W）
w=
4.0000+0.0000i1.5000+0.0000i1.0000+0.0000i1.0000+0.0000i
1.0000+0.0000i1.5000+0.0000i%虚部都是0，说明是实数；
w（n）={0,0,4,3,2,1}
图3.5
升级会员 