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

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

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

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

学生姓名：闫春遐

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

年级专业：2008级自动化系

学 号：00824049

指导教师：王亮

完成日期：2010年9月27日
实验二 快速傅里叶变换（FFT）及其应用
一、实验目的
（1）在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。
（2）应用FFT对典型信号进行频谱分析。
（3）了解应用FFT进行信号频谱分析过程可能出现的问题，以便在实际中正确应用FFT。
（4）应用FFT实现序列的线性卷积和相关。
二、实验内容
实验中用到的信号序列：
a） 高斯序列

b） 衰减正弦序列

c） 三角波序列

d） 反三角波序列

上机实验内容：
（1）观察高斯序列的时域和幅频特性，固定信号中参数，改变的值，使分别等于2、4、8，观察他们的时域和幅频特性，了解当取不同值时，对信号的时域和幅频特性的影响；固定，改变，使分别等于8、13、14，观察参数变化对信号序列的时域及幅频特性的影响，注意等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。
解答:
>> 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）观察衰减正弦序列的时域和幅频特性，，，检查普峰出现的位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变，使分别等于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）观察三角波和反三角波的时域和幅频特性，用点FFT分析信号序列和的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。
在和末尾补零，用点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)');

时，和的幅频特性相同，在和末尾补零，用点FFT分析这两个信号的幅频特性时，它们还有相同之处，即当取4的整数倍时对应幅值相等。
分析：
点FFT分析信号的幅频特性：
点FFT分析信号的幅频特性：

由上两式可知，当k2=4k1时，两个信号的对应频率幅值相等，即对信号末尾补零加长整数个周期可以对原信号达到细化频谱的作用。
（4）一个连续时间信号含两个频率分量，经采样得

已知，分别为1/16和1/64，观察其频谱；当时，不变，其结果有何不同，为什么？
解答：
>> 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');

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

当不等于整数时，则信号频谱会发生泄漏。
（5）用FFT分别计算（）和（）的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点的随机序列，并用和做线性卷积，观察卷积前后频谱的变化。要求将分成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;
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;
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分别计算（）和（）的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分别计算（）和（）的自相关函数。
解答：
>> 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('幅度');

第18页 共19页