实验二-快速傅里叶变换(FFT)及其应用

实验二-快速傅里叶变换(FFT)及其应用

《数字信号处理》课程 （2010-2011学年第1学期 ）

成绩：

实验二 快速傅里叶变换（FFT）及其应用

学生姓名：闫春遐

所在院系：电子信息工程学院自动化系

年级专业：2008级自动化系

学 号：00824049

指导教师：王亮

完成日期：2010年9月27日

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

实验二 快速傅里叶变换（FFT）及其应用

一、实验目的

（1）在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。

（2）应用FFT对典型信号进行频谱分析。

（3）了解应用FFT进行信号频谱分析过程可能出现的问题，以便在实际中正确应用FFT。

（4）应用FFT实现序列的线性卷积和相关。

二、实验内容

实验中用到的信号序列： a）

高斯序列

2(np)qxa(n)e00n15

其他b） 衰减正弦序列

eansin(2fn)0n15 xb(n)

0其他c）

三角波序列

0n3nxc(n)8n4n70其他d） 反三角波序列

4n0n3xd(n)n44n70其他上机实验内容：

（1）观察高斯序列的时域和幅频特性，固定信号xa(n)中参数p8，改变q的值，使q分别等于2、4、8，观察他们的时域和幅频特性，了解当q取不同值时，对信号的时域和幅频特性的影响；固定q8，改变p，使p分别等于8、13、

第1页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

14，观察参数p变化对信号序列的时域及幅频特性的影响，注意p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。

解答: >> n=0:1:15;

>> xn=exp(-(n-8).^2/2);

>> subplot(1,2,1);stem(n,xn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

>> xn=exp(-(n-8).^2/4);

>> subplot(1,2,1);stem(n,xn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

第2页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

>> xn=exp(-(n-8).^2/8);

>> subplot(1,2,1);stem(n,xn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

>> xn=exp(-(n-13).^2/8);

>> subplot(1,2,1);stem(n,xn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

>> xn=exp(-(n-14).^2/8);

>> subplot(1,2,1);stem(n,xn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

第3页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

随着q值的增大，时域信号幅值变化缓慢，频域信号频谱泄露程度减小。 随着p的增大，时域信号幅值不变，会在时间轴移位。

（2）观察衰减正弦序列xb(n)的时域和幅频特性，a0.1，f0.0625，检查普峰出现的位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使

f分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和普峰出现的位

置，有无混叠和泄漏现象？说明产生现象的原因。

解答： >> n=0:1:15;

>> xn=exp(-0.1*n).*sin(2*pi*0.0625*n);

>> subplot(1,2,1);stem(n,xn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

第4页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

>> xn=exp(-0.1*n).*sin(2*pi*0.4375*n);

>> subplot(1,2,1);stem(n,xn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

>> xn=exp(-0.1*n).*sin(2*pi*0.5625*n);

>> subplot(1,2,1);stem(n,xn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

（3）观察三角波和反三角波的时域和幅频特性，用N8点FFT分析信号序列xc(n)和xd(n)的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。

在xc(n)和xd(n)末尾补零，用N32点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？两种情况下的FFT频谱还有相同之处吗？这些变

第5页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

化说明了什么？

解答：

>> for n=0:1:3 xcn(n+1)=n; end;

>> for n=4:1:7 xcn(n+1)=8-n; end; >> xcn xcn =

0 1 2 3 4 3 2 1

>> n=0:1:7;

>> subplot(1,2,1);stem(n,xcn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xcn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

>> for n=0:1:3 xdn(n+1)=4-n; end;

第6页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

>> for n=4:1:7 xdn(n+1)=n-4; end; >> xdn xdn =

4 3 2 1 0 1 2 3

>> n=0:1:7;

>> subplot(1,2,1);stem(n,xdn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xdn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

>> xcn=[xcn,zeros(1,24)]; >> n=0:1:31;

>> subplot(1,2,1);stem(n,xcn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xcn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

第7页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

>> xdn=[xdn,zeros(1,24)]; >> n=0:1:31;

>> subplot(1,2,1);stem(n,xdn);xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xdn);xk1=abs(xk1);

>> subplot(1,2,2);stem(n,xk1);xlabel('k');ylabel('X(k)');

N8时，xc(n)和xd(n)的幅频特性相同，在xc(n)和xd(n)末尾补零，用

N32点FFT分析这两个信号的幅频特性时，它们还有相同之处，即当k取4的整数倍时对应幅值相等。

分析：

N1n02)nk1NN8点FFT分析信号的幅频特性：XN(k1)x(n)*eN32点FFT分析信号的幅频特性：

4N1j(

X4N(k2)x(n)*en0j(2)nk24Nx(n)*en0N1j(2)nk24N

第8页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

由上两式可知，当k2=4k1时，两个信号的对应频率幅值相等，即对信号末尾补零加长整数个周期可以对原信号达到细化频谱的作用。

（4）一个连续时间信号含两个频率分量，经采样得

x(n)sin[20.125n]cos[2(0.125f)n]n0,1,,N1

已知N16，f分别为1/16和1/64，观察其频谱；当N128时，f不变，其结果有何不同，为什么？

解答： >> n=0:1:15;

>> x1n=sin(2*pi*0.125*n)+cos(2*pi*(0.125+1/16)*n); >> xk1=fft(x1n);xk1=abs(xk1);

>>subplot(1,2,1);stem(n,xk1);xlabel('k');ylabel('X(k)');legend('f=1/16');

>> x2n=sin(2*pi*0.125*n)+cos(2*pi*(0.125+1/64)*n); >> xk2=fft(x2n);xk2=abs(xk2);

>>subplot(1,2,2);stem(n,xk2);xlabel('k');ylabel('X(k)');legend('f=1/64');

>> n=0:1:127;

>> x1n=sin(2*pi*0.125*n)+cos(2*pi*(0.125+1/16)*n); >> xk1=fft(x1n);xk1=abs(xk1);

>> stem(n,xk1);xlabel('k');ylabel('X(k)');legend('f=1/16');

第9页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

>> x2n=sin(2*pi*0.125*n)+cos(2*pi*(0.125+1/64)*n); >> xk2=fft(x2n);xk2=abs(xk2);

>> stem(n,xk2);xlabel('k');ylabel('X(k)');legend('f=1/64');

分析：

由于离散傅里叶变换的选频性质：

x(n)ejqwonwo2/N

qkqk

N1ej2(qk)X(k)DFT[x(n)]1ej2(qk)/N0当q不等于整数时，则信号频谱会发生泄漏。

（5）用FFT分别计算xa(n)（p8,q2）和xb(n)（a0.1,f0.0625）的16点循环卷积和线性卷积。

解答： >> n=0:1:15;

>> xan=exp(-(n-8).^2/2);

>> xbn=exp(-0.1*n).*sin(2*pi*0.0625*n);

第10页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

>> subplot(4,1,1);stem(n,xan);xlabel('n');ylabel('xa(n)'); >> subplot(4,1,2);stem(n,xbn);xlabel('n');ylabel('xb(n)'); >> xak=fft(xan);xbk=fft(xbn);x1k=xak.*xbk; >> x1n=ifft(x1k);

>>subplot(4,1,3);stem(n,x1n);xlabel('n');ylabel('x1(n)');legend('循环卷积');

>> x2n=conv(xan,xbn); >> m=0:1:length(x2n)-1;

>>subplot(4,1,4);stem(m,x2n);xlabel('n');ylabel('x2(n)');legend('线性卷积');

（6）产生一512点的随机序列xe(n)，并用xc(n)和xe(n)做线性卷积，观察卷积前后xe(n)频谱的变化。要求将xe(n)分成8段，分别采用重叠相加法和重叠保留法。

解答：

第11页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

在编辑调试窗中编写程序： function yy=xeni(N2,xen,i) for n=N2*i:1:N2*(i+1)-1 xeni(n-N2*i+1)=xen(n+1); end yy=xeni;

将上述文件存盘，文件名为xeni.m。

function yy=xenni(N1,N2,xen,i) for n=N2*i:1:N1+N2*(i+1)-2 xeni(n-N2*i+1)=xen(n+1); end yy=xeni;

将上述文件存盘，文件名为xenni.m。

function t=shiftmm(a,n) m=length(n); for i=1:1:a;

for j=m+i-1:-1:1 n(j+1)=n(j); end; end; for i=1:1:a n(i)=0; end; t=n;

将上述文件存盘，文件名为shiftmm.m。 退回到指令窗：

>> xcn=[0 1 2 3 4 3 2 1];xen=rand(1,512);

第12页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

>> qqqqq=conv(xcn,xen);

>> stem([0:1:518],qqqqq);xlabel('n');ylabel('幅度');

>> N1=length(xcn);N2=length(xen)/8; >> xcn=[xcn zeros(1,N2-1)]; >> xck=fft(xcn); >> for i=1:1:8

xenii=xeni(N2,xen,i-1); xenii=[xenii zeros(1,N1-1)]; xeki=fft(xenii); yki=xck.*xeki; yni=ifft(yki); y(i,:)=yni; end;

>> for i=0:1:7 for j=0:1:i*N2-1

ynii(i+1,[0+1:1:i*N2-1+1])=0; end;

for j=i*N2:1:N1+(i+1)*N2-2

ynii(i+1,[i*N2+1:1:N1+(i+1)*N2-2+1])=y(i+1,:); end;

for j=N1+(i+1)*N2-1:1:N1+8*N2-2

ynii(i+1,[N1+(i+1)*N2-1+1:1:N1+8*N2-2+1])=0; end;

第13页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

end;

>> yn=zeros(1,N1+8*N2-1); >> for i=1:1:8 yn=yn+ynii(i,:); end;

>> n=0:1:N1+8*N2-2;

>> stem(n,yn);xlabel('n');ylabel('幅度');legend('重叠相加法');

>> xen21=shiftmm(N1-1,xen); >> for i=1:1:8

xen2i(i,:)=xenni(N1,N2,xen21,i-1); end;

>> for i=1:1:8

xek2i=fft(xen2i(i,:)); yk2i=xck.*xek2i; yn2i=ifft(yk2i); y2(i,:)=yn2i;; end;

>> y2(:,1:N1-1)=[;;;;;;;;]; >> n2=0:1:8*N2-1;

>> stem(n2,[y2(1,:) y2(2,:) y2(3,:) y2(4,:) y2(5,:) y2(6,:) y2(7,:) y2(8,:)]);xlabel('n');ylabel('幅度');legend('重叠保留法');

第14页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

（7）用FFT分别计算xa(n)（p8,q2）和xb(n)（a0.1,f0.0625）的16点循环相关和线性相关，问一共有多少种结果，它们之间有何异同点。

解答： 1）求线性相关 >> n=0:1:15;

>> xan=exp(-(n-8).^2/2);

>> xbn=exp(-0.1*n).*sin(2*pi*0.0625*n); >> k=length(xbn);

>> xan1=[xan zeros(1,k-1)]; >> xbn1=[xbn zeros(1,k-1)]; >> xak=fft(xan1); >> xbk=fft(xbn1);

>> rm=real(ifft(conj(xak).*xbk)); >> rm1=[rm(k+1:2*k-1) rm(1:k)]; >> m=(-k+1):(k-1);

>> stem(m,rm1);xlabel('n');ylabel('幅度');legend('线性相关');

第15页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

2）求循环相关 >> n=0:1:15;

>> xan=exp(-(n-8).^2/2);

>> xbn=exp(-0.1*n).*sin(2*pi*0.0625*n); >> k=length(xbn); >> xak=fft(xan); >> xbk=fft(xbn);

>> rm=real(ifft(conj(xak).*xbk));

>> stem(n,rm);xlabel('n');ylabel('幅度');legend('循环相关');

第16页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

（8）用FFT分别计算xa(n)（p8,q2）和xb(n)（a0.1,f0.0625）的自相关函数。

解答： >> n=0:1:15;

>> xan=exp(-(n-8).^2/2); >> k=length(xan); >> xak=fft(xan,2*k);

>> rm=real(ifft(conj(xak).*xak)); >> rm=[rm(k+2:2*k) rm(1:k)]; >> m=(-k+1):(k-1);

>> stem(m,rm);xlabel('m');ylabel('幅度');

（2） >> n=0:1:15;

>> xbn=exp(-0.1*n).*sin(2*pi*0.0625*n); >> k=length(xbn); >> xbk=fft(xbn,2*k);

>> rm=real(ifft(conj(xbk).*xbk)); >> rm=[rm(k+2:2*k) rm(1:k)]; >> m=(-k+1):(k-1);

第17页 共20页

实验二 快速傅里叶变换（FFT）及其应用 闫春遐 00824049

>> stem(m,rm);xlabel('m');ylabel('幅度');

第18页 共20页