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

《数字信号处理》课程

（2010-2011学年第1学期）成绩：

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

学生姓名：闫春遐

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

年级专业：2008级自动化系

学号：00824049

指导教师：王亮

完成日期：2010年9月27日

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

一、实验目的

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

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

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

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

二、实验内容

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

()015()0

n p q

a e

n x n --??≤≤=???其他

b ）衰减正弦序列

sin(2)015

()0an b e fn n x n π-?≤≤=??

其他

c ）三角波序列

03()847

0c n

n x n n n ≤≤??

=-≤≤???

其他

d ）反三角波序列

403()447

0d n n x n n n -≤≤??

=-≤≤???

其他

上机实验内容：

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

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)');

>> 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)');

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

（2）观察衰减正弦序列()b x n 的时域和幅频特性，0.1a =，0.0625f =，检查普峰出现的位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变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)');

>> 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）观察三角波和反三角波的时域和幅频特性，用8N =点FFT 分析信号序列()c x n 和()d x n 的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。

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

化说明了什么？

解答：

>> 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;

>> 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)');

>> 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)');

8N =时，()c x n 和()d x n 的幅频特性相同，在()c x n 和()d x n 末尾补零，用

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

分析：

8N =点FFT 分析信号的幅频特性：1121(

()()*N N j nk N

n X k x n e π

-==∑

32N =点FFT 分析信号的幅频特性：

22422241

(

)(

)440

()()*()*N N N j nk j nk N

n n X k x n e

x n e

--===

=∑∑

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

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

()sin[20.125]cos[2(0.125)]0,1,,1x n n f n n N ππ=++?=???-g g

已知16N =，f ?分别为1/16和1/64，观察其频谱；当128N =时，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');

>> 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');

分析：

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

()2/o jqw n

o x n e w N

π==

2()

2()/1()[()]0

1j q k j q k N N q k e X k DFT x n q k

e ππ--=?-===?

≠-?

当q 不等于整数时，则信号频谱会发生泄漏。

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

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

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

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

>> 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点的随机序列()e x n ，并用()c x n 和()e x n 做线性卷积，观察卷积前后()e x n 频谱的变化。要求将()e x n 分成8段，分别采用重叠相加法和重叠保留法。

解答：

在编辑调试窗中编写程序：

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;

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);

>> 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;

>> 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('重叠保留法');

（7）用FFT 分别计算()a x n （8,2p q ==）和()b x n （0.1,0.0625a f ==）的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('线性相关');

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('循环相关');

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

解答： >> 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);

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