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

第1页/共17页实验报告课程名称：信号分析与处理指导老师成绩：实验名称：离散傅里叶变换和快速傅里叶变换实验类型：同组学生姓名：一、实验目的和要求（必填）二、实验内容和原理（必填）三、主要仪器设备（必填）四、操作方法和实验步骤五、实验数据记录和处理六、实验结果与分析（必填）七、讨论、心得一、实验目的和要求1.掌握DFT的原理和实现2.掌握FFT的原理和实现，掌握用FFT对连续信号和离散信号进行谱分析的方法。二、实验内容和原理2.1DTFT和DFT序列x(n)的离散事件傅里叶变换（DTFT）表示为：njnjenxeX)()(，如果x(n)为因果有限长序列，n=0,1,...,N-1，则x(n)的DTFT表示为：njNnjenxeX10)()(，x(n)的离散傅里叶变换（DFT）表达式为：)1,...,1,0()()(210NkenxkXnkNjNn，序列的N点DFT是DTFT在[0,2π]上的N点等间隔采样，采样间隔为2π/N。通过DFT，可以完成由一组有限个信号采样值x(n)直接计算得到一组有限个频谱采样值X(k)。X(k)的幅度谱为)()()(22kXkXkXIR，XR(k)和XI(k)分别为X(k)的实部和虚部。X(k)的相位谱为)()(arctan)(kXkXkRI。第2页/共17页离散傅里叶反变换（IDFT）定义为)1,...,1,0()(1)(210NnekXNnxnkNjNn。2.2FFT快速傅里叶变换（FFT）是DFT的快速算法，它减少了DFT的运算量，使数字信号的处理速度大大提高。三、主要仪器设备PC一台，matlab软件四、实验内容4.1第一题求有限长离散时间信号x(n)的离散时间．．傅里叶变换（DTFT）X(ejΩ)并绘图。（1）已知其他0221)(nnx；（2）已知1002)(nnxn。4.1.1理论分析1)由DTFT计算式，252.52.520.50.52e1eeesin(2.5)()()ee1eeesin(0.5)jjjjjnjnjjjnnXxnX（Ω）是实数，可以直接作出它的图像。第3页/共17页Figure1X（Ω）曲线2)由DTFT计算式：1111111110012e2e()()e2e12e12ejjjnnjnjjnnXxn112()12ejX可以发现X（Ω）周期为2π；而X（Ω）的相位在2π周期内有约十次振荡。4.1.2编程计算作图编写一个计算DTFT的函数。functionDTFT(x,n1,n2)w=-2*pi:2*pi/1000:2*pi;%表示ΩX=zeros(size(w));fori=n1:n2%DTFT计算式X=X+x(i-n1+1)*exp((-1)*j*w*i);endangle(X);subplot(2,1,1);plot(w,abs(X),'r');xlabel('\Omega');ylabel('|X(\Omega)|');holdon;%作幅频图subplot(2,1,2);plot(w,angle(X),'b');xlabel('\Omega');ylabel('\angle(\Omega)');%作相频图end输入序列x，和n的取值范围，即可计算其DTFT。第4页/共17页1)输入：x=[11111];DTFT(x,-2,2);（因为X（Ω）是实数，所以实际计算过程中对相频曲线取了绝对值）结果：-8-6-4-202468-1012345|X()|-8-6-4-202468-10123angle()Figure2X（Ω）的频谱可以看出，X（Ω）的相位只有0和π两种取值，X（Ω）是实函数，而且其幅度频谱与理论计算得到的相同。2)输入：n=0:10;x=2.^n;DTFT(x,0,10)结果：-15-10-505105001000150020002500|X()|-15-10-50510-4-2024angle()Figure3第1题（2）中X（Ω）的频谱第5页/共17页4.2第二题已知有限长序列x(n)={8,7,9,5,1,7,9,5},试分别采用DFT和FFT求其离散傅里叶变换X(k)的幅度、相位图。4.2.1理论分析由FFT蝶形运算得到，X(k)=51,7,-9-j4,7,3,7,-9+j4,74.2.2编程计算作图1.DFT编写一个计算DTFT的函数。DFT（序列x，长度N）functionDFT(x,N)k=0:N-1;X=zeros(size(k));forn=0:N-1X=X+x(n+1)*exp((-1)*j*2*pi/N*n*k);%DFT计算式endsubplot(2,1,1);stem(k,abs(X),'.');xlabel('k');ylabel('|X(k)|');holdon;%幅频图subplot(2,1,2);stem(k,angle(X),'*');xlabel('k');ylabel('Angle(k)');%相频图end输入：x=[87951795];DFT(x,8);结果：-10123456780204060k|X(k)|-1012345678-505kAngle(k)Figure4第2题DFT结果第6页/共17页2.FFT编写一个利用matlab自带函数计算FFT并绘图的函数FFT1（序列x，长度N）functionFFT1(x,N)X=fft(x,8);%用自带的fft函数计算k=0:N-1;subplot(2,1,1);stem(k,abs(X),'.');xlabel('k');ylabel('|X(k)|');holdon;%幅频图subplot(2,1,2);stem(k,angle(X),'*');xlabel('k');ylabel('Angle(k)');%相频图end输入：x=[87951795];FFT1(x,8);结果：-10123456780204060k|X(k)|-1012345678-505kAngle(k)Figure5第2题FFT结果因为FFT只是DFT的一种快速算法，所以FFT的结果与DFT结果相同。DFT和FFT的结果，符合理论计算得到的，X(k)=51,7,-9-j4,7,3,7,-9+j4,7。4.3第三题已知连续时间信号x(t)=3cos8πt,X(ω)=)]8()8([3，该信号从t=0开始以采样周期Ts=0.1s进行采样得到序列x(n)，试选择合适的采样点数，分别采用DFT和FFT求其离散傅里叶变换X(k)的幅度、相位图，并将结果与X(k)的幅度、相位图，并将结果与X(ω)相比较。4.3.1理论分析1.原信号的频谱：X(ω)=)]8()8([3，只在±8π不为0.且在±8π处相位为0。第7页/共17页2.采样角频率Ωs=20π＞2×8π，满足采样定理。3.采样后的信号，为X(ω)以20π为周期的延拓。所以只在（±8+20k）π（k为任意整数）处不为0.如取区间[0,20π]内，只有8π和16π处不为0。进行N点DFT后，将20π的区间映射为[0,N]区间。理想情况下仅在n=820N和1220N两处不为0。4.(n)cos(8nTs)cos(0.8)xn，周期为5，所以取采样点数为5的倍数时，不会发生泄漏；而采样点数不是5的倍数时，则会发生泄漏。-40-30-20-10010203040012345678910X()33-88Figure6原始信号的频谱X（Ω）4.3.2编程计算作图编写一个获得信号的N点样本的函数sample（点数N）functionx=sample(N)t=0:0.1:(N-1)*0.1;%0.1s为间隔x=3*cos(8*pi*t);%x即采样结果。End输入:X=sample(N);FFT1(X,N);即可获得采样N点的频谱图。因为DFT结果与FFT是完全一样的，所以这里只使用FFT作图。取采样点数N=51620104获得以下频谱图：第8页/共17页00.511.522.533.5402468k|X(k)|00.511.522.533.54-2024kAngle(k)Figure7N=50510150102030k|X(k)|051015-2-1012kAngle(k)Figure8N=1602468101214161820010203040k|X(k)|02468101214161820-4-2024kAngle(k)Figure9N=20第9页/共17页0102030405060020406080k|X(k)|0102030405060-2-1012kAngle(k)Figure10N=104可以看出，N=5和20时，由于是周期的整数倍，频谱只有两条谱线，且满足前面理论计算得出的公式n=820N和1220N，没有发生泄漏，且这两条谱线对应的相位是0.所以频谱与原信号频谱在形式上时相同的。而N=16和N=54时，则都发生了频谱泄漏，频谱与原信号频谱就很不同了。但相比之下N=54时谱线更加接近原谱线。验证了“为减小泄漏误差，如果待分析的信号实现不知道确切周期，则截取较长时间长度的样点进行分析”这个说法。同时也可以发现，虽然幅频图中显示幅值为0，但相频图中相应的位置仍有谱线。这可能是matlab浮点运算造成的误差，即本来为0处其实是一个非常小的复数，所以仍有一定相位。4.4第四题4.4.1理论分析若噪声信号较小，则采样后的频谱仍能较准确地反映原信号的特征。4.4.2编程计算对原采样程序稍加改编，加入一个噪声信号p*randn(1,N)。p表示噪声信号的强度。functionx=samplenoise(N,p)t=0:0.1:(N-1)*0.1;第10页/共17页x=3*cos(8*pi*t)+p*randn(1,N);end取采样点数N=20进行分析。输入：X=samplenoise(20,p);%取P=1和10两种情况。FFT1(X,20);024681012141618200102030k|X(k)|02468101214161820-4-2024kAngle(k)Figure11N=20噪声较小（p=1）02468101214161820050100k|X(k)|02468101214161820-4-2024kAngle(k)Figure12N=20噪声较大（p=10）可见，较小的噪声对信号的频谱的影响不大，仍能较精确地获得频谱图。而噪声较大时则很难准确获得原信号的频谱。4.5第五题3.5已知序列)3()2(2)1(3)(4)(nnnnnx，X(k)是x(n)的6点DFT，设kNjkNeW2。（1）若有限长序列y(n)的6点DFT是)()(46kXWkYk，求y(n)。第11页/共17页（2）若有限长序列w(n)的6点DFTW(k)是)(kX的实部，求w(n)。（3）若有限长序列q(n)的3点DFT是)2()(kXkQ，k=0,1,2，求q(n)。由题意得到：x(n)=4,3,2,1,0,0n=(0,1,2,3,4,5)4.5.1理论分析1)由DFT的性质可以得到，如果y（n）的DFT为46()()kYkWXk，那么y(n)就是x(n)圆周左移4位得到的。所以y（n）=0,0,4,3,2,12)由题意，X(k)=103.5-j4.33012.5-8.66j22.5+j0.8663.5+j4.3301取实部，则W(k)=103.52.52.523.5按照IDFT计算式)1,...,1,0()(1)(210NnekXNnxnkNjNn计算得到w(n)=41.5111.5由DFT的性质，因为W(k)是实数，所以对应的w(n)也是实数。3)由题意，Q(k)=102