数字信号处理吴镇杨实验三.docx
- 文档编号:4139975
- 上传时间:2022-11-28
- 格式:DOCX
- 页数:17
- 大小:263.54KB
数字信号处理吴镇杨实验三.docx
《数字信号处理吴镇杨实验三.docx》由会员分享,可在线阅读,更多相关《数字信号处理吴镇杨实验三.docx(17页珍藏版)》请在冰豆网上搜索。
数字信号处理吴镇杨实验三
实验三快速傅里叶变换及其应用
04011344王晨
一、实验目的
(1)在理论学习的基础上,通过本实验,加深对FFT的理解,熟悉MATLAB中的有关函数。
(2)应用FFT对典型信号进行频谱分析的方法。
(3)了解应用FFT进行信号频谱分析过程中可能出现的问题,以便在实际中正确应用FFT。
(4)应用FFT实现两个序列的线性卷积和相关
二、实验原理
在各种信号序列中,有限长序列信号处理占有很重要地位,对有限长序列,我们可以使用离散Fourier变换(DFT)。
这一变换不但可以很好的反映序列的频谱特性,而且易于用快速算法在计算机上实现,当序列x(n)的长度为N时,它的DFT定义为:
反变换为:
有限长序列的DFT是其Z变换在单位圆上的等距采样,或者说是序列Fourier变换的等距采样,因此可以用于序列的谱分析。
FFT并不是与DFT不同的另一种变换,而是为了减少DFT运算次数的一种快速算法。
它是对变换式进行一次次分解,使其成为若干小点数的组合,从而减少运算量。
常用的FFT是以2为基数的,其长度
。
它的效率高,程序简单,使用非常方便,当要变换的序列长度不等于2的整数次方时,为了使用以2为基数的FFT,可以用末位补零的方法,使其长度延长至2的整数次方。
(一)、在运用DFT进行频谱分析的过程中可能产生三种误差:
(1) 混叠
序列的频谱时被采样信号的周期延拓,当采样速率不满足Nyquist定理时,就会发生频谱混叠,使得采样后的信号序列频谱不能真实的反映原信号的频谱。
避免混叠现象的唯一方法是保证采样速率足够高,使频谱混叠现象不致出现,即在确定采样频率之前,必须对频谱的性质有所了解,在一般情况下,为了保证高于折叠频率的分量不会出现,在采样前,先用低通模拟滤波器对信号进行滤波。
(2) 泄漏
实际中我们往往用截短的序列来近似很长的甚至是无限长的序列,这样可以使用较短的DFT来对信号进行频谱分析,这种截短等价于给原信号序列乘以一个矩形窗函数,也相当于在频域将信号的频谱和矩形窗函数的频谱卷积,所得的频谱是原序列频谱的扩展。
泄漏不能与混叠完全分开,因为泄漏导致频谱的扩展,从而造成混叠。
为了减少泄漏的影响,可以选择适当的窗函数使频谱的扩散减至最小。
(3) 栅栏效应
DFT是对单位圆上Z变换的均匀采样,所以它不可能将频谱视为一个连续函数,就一定意义上看,用DFT来观察频谱就好像通过一个栅栏来观看一个图景一样,只能在离散点上看到真实的频谱,这样就有可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所拦住,不能别我们观察到。
减小栅栏效应的一个方法就是借助于在原序列的末端填补一些零值,从而变动DFT的点数,这一方法实际上是人为地改变了对真实频谱采样的点数和位置,相当于搬动了每一根“尖桩栅栏”的位置,从而使得频谱的峰点或谷点暴露出来。
(二)、用FFT计算线性卷积
用FFT可以实现两个序列的圆周卷积。
在一定的条件下,可以使圆周卷积等于线性卷积。
一般情况,设两个序列的长度分别为N1和N2,要使圆周卷积等于线性卷积的充要条件是FFT的长度:
N≥N1+N2
对于长度不足N的两个序列,分别将他们补零延长到N。
当两个序列中有一个序列比较长的时候,我们可以采用分段卷积的方法。
有两种方法:
a.重叠相加法。
将长序列分成与短序列相仿的片段,分别用FFT对它们作线性卷积,再将分段卷积各段重叠的部分相加构成总的卷积输出。
b.重叠保留法。
这种方法在长序列分段时,段与段之间保留有互相重叠的部分,在构成总的卷积输出时只需将各段线性卷积部分直接连接起来,省掉了输出段的直接相加。
(三)、用周期图法(平滑周期图的平均法)对随机信号作谱分析
实际中许多信号往往既不具有有限能量,由非周期性的。
无限能量信号的基本概念是随机过程,也就是说无限能量信号是一随机信号。
周期图法是随机信号作谱分析的一种方法,它特别适用于用FFT直接计算功率谱的估值。
将长度为N的实平稳随机序列的样本x(n)再次分割成K段,每段长度为L,即L=N/K。
每段序列仍可表示为:
xi(n)=x(n+(i-1)L),0≤n≤L-1,1≤i≤K
但是这里在计算周期图之前,先用窗函数w(n)给每段序列xi(n)加权,K个修正的周期图定义为
其中U表示窗口序列的能量,它等于
在此情况下,功率谱估计量可表示为
三、实验程序代码及结果讨论
实验中用到的信号序列:
a) Gaussian序列
b) 衰减正弦序列
c) 三角波序列
d) 反三角波序列
上机实验内容:
(1)、观察高斯序列的时域和幅频特性,固定信号xa(n)中参数p=8,改变q的值,使q分别等于2,4,8,观察它们的时域和幅频特性,了解当q取不同值时,对信号序列的时域幅频特性的影响;固定q=8,改变p,使p分别等于8,13,14,观察参数p变化对信号序列的时域及幅频特性的影响,观察p等于多少时,会发生明显的泄漏现象,混叠是否也随之出现?
记录实验中观察到的现象,绘出相应的时域序列和幅频特性曲线。
高斯序列定义:
function[Xa,Fa]=gauss(p,q)
n=[0:
15];
Xa(n+1)=exp(-(n+1-p).^2./q);
F=fft(Xa);
Fa=abs(F);
end
clearall;
[Xa1,Fa1]=gauss(8,2);
k=0:
15;
subplot(5,2,1);stem(k,Xa1);hold;plot(k,Xa1);
xlabel('n');ylabel('时域特性');text(10,0.5,'p=8,q=2');
subplot(5,2,2);stem(k,Fa1);hold;plot(k,Fa1);
xlabel('n');ylabel('幅频特性');text(8,3,'p=8,q=2');
[Xa2,Fa2]=gauss(8,4);
subplot(5,2,3);stem(k,Xa2);hold;plot(k,Xa2);
xlabel('n');ylabel('时域特性');text(10,0.5,'p=8,q=4');
subplot(5,2,4);stem(k,Fa2);hold;plot(k,Fa2);
xlabel('n');ylabel('幅频特性');text(8,3,'p=8,q=4');
[Xa3,Fa3]=gauss(8,8);
subplot(5,2,5);stem(k,Xa3);hold;plot(k,Xa3);
xlabel('n');ylabel('时域特性');text(10,0.5,'p=8,q=8');
subplot(5,2,6);stem(k,Fa3);hold;plot(k,Fa3);
xlabel('n');ylabel('幅频特性');text(8,3,'p=8,q=8');
[Xa4,Fa4]=gauss(13,8);
subplot(5,2,7);stem(k,Xa4);hold;plot(k,Xa4);
xlabel('n');ylabel('时域特性');text(10,0.5,'p=13,q=8');
subplot(5,2,8);stem(k,Fa4);hold;plot(k,Fa4);
xlabel('n');ylabel('幅频特性');text(8,3,'p=13,q=8');
[Xa5,Fa5]=gauss(14,8);
subplot(5,2,9);stem(k,Xa5);hold;plot(k,Xa5);
xlabel('n');ylabel('时域特性');text(10,0.5,'p=14,q=8');
subplot(5,2,10);stem(k,Fa5);hold;plot(k,Fa5);
xlabel('n');ylabel('幅频特性');text(8,3,'p=14,q=8');
实验现象分析:
由高斯序列表达式知n=p为其对称轴;
当p取固定值时,时域序列图以n=p为其对称轴,当q取值越来越大时,波形变化更加缓和,即向对称轴两边下降越加平缓。
q值越大,高频分量越少,混叠效应减弱。
由于原序列相当于是对称截取,且窗口包含序列主要部分,故泄露现象并不明显。
当q值固定不变,p变化时,时域对称轴右移,截取窗口不能捕捉信号主要部分,开始无法代表一个周期,泄漏现象也越来越明显,泄露也产生一定的频谱混叠。
p=14时的泄漏现象最为明显,混叠也相应产生;
(2)、观察衰减正弦序列xb(n)的时域和幅频特性,a=0.1,f=0.0625,检查谱峰出现位置是否正确,注意频谱的形状,绘出幅频特性曲线,改变f,使f分别等于0.4375和0.5625,观察这两种情况下,频谱的形状和谱峰出现位置,有无混叠和泄漏现象?
说明产生现象的原因。
衰减正弦序列定义:
function[Xb,Fb]=downsin(a,f)
n=[0:
15];
Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n);
F=fft(Xb);
Fb=abs(F);
clearall;
k=0:
15;
[Xb,Fb]=downsin(0.1,0.0625);
subplot(3,2,1);stem(k,Xb);hold;plot(k,Xb);
xlabel('n');ylabel('时域特性');text(8,0.5,'a=0.1,f=0.0625');
subplot(3,2,2);stem(k,Fb);hold;plot(k,Fb);
xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.0625');
[Xb1,Fb1]=downsin(0.1,0.4375);
subplot(3,2,3);stem(k,Xb1);hold;plot(k,Xb1);
xlabel('n');ylabel('时域特性');text(8,0.5,'a=0.1,f=0.4375');
subplot(3,2,4);stem(k,Fb1);hold;plot(k,Fb1);
xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.4375');
[Xb2,Fb2]=downsin(0.1,0.5625);
subplot(3,2,5);stem(k,Xb2);hold;plot(k,Xb2);
xlabel('n');ylabel('时域特性');text(8,0.5,'a=0.1,f=0.5625');
subplot(3,2,6);stem(k,Fb2);hold;plot(k,Fb2);
xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.5625');
实验现象分析:
当f=f1/fs=1/16=0.0625时,谱峰位置出现正确。
时间长度为16Ts,即采样频率为16f1,满足采样定理。
由于截取长度为周期的整数倍(一倍),此处不引入频谱泄漏。
但是由于截断效应,相当于加一矩形窗,其频谱中旁瓣导致频率分量扩散,因此有一定程度的频谱泄漏。
增加的部分高频分量会造成低频成分的混叠效应,但是由截断带来的高频旁瓣较小,故混叠效应较弱。
当f=f2/fs=7/16=0.4375时,fs=(16/7)f2,仍满足采样定理,谱峰位置正确,且为整数周期采样。
同上,较弱的频谱泄漏及混叠效应来自于截断效应。
当f=f2/fs=9/16=0.5625时,fs=(16/9)f2,不满足采样定理,故产生混叠,此时谱峰位置不正确。
混叠后的频谱由图知与第二种情况频谱相同,原因是超过1/2fs的高频分量叠加到低频。
此时也产生了一定的频谱泄漏。
(3)、观察三角波和反三角波序列的时域和幅频特性,用N=8点FFT分析信号序列xc(n)和xd(n)的幅频特性,观察两者的序列形状和频谱曲线有什么异同?
绘出两序列及其幅频特性曲线。
在xc(n)和xd(n)末尾补零,用N=32点FFT分析这两个信号的幅频特性,观察幅频特性发生了什么变化?
两情况的FFT频谱还有相同之处吗?
这些变化说明了什么?
程序代码如下:
clearall;
n=[0:
3];k=[1:
8];
Xc(n+1)=n;Xc(n+5)=4-n;
Xd(n+1)=4-n;Xd(n+5)=n;
subplot(2,2,1);stem(k-1,Xc);hold;plot(k-1,Xc);
xlabel('n');ylabel('时域特性');text(1,3,'三角波');
subplot(2,2,2);stem(k-1,abs(fft(Xc)));hold;plot(k-1,abs(fft(Xc)));
xlabel('k');ylabel('幅频特性');text(4,10,'三角波');
subplot(2,2,3);stem(k-1,Xd);hold;plot(k-1,Xd);
xlabel('n');ylabel('时域特性');text(3,3,'反三角波');
subplot(2,2,4);stem(k-1,abs(fft(Xd)));hold;plot(k-1,abs(fft(Xd)));
xlabel('k');ylabel('幅频特性');text(4,10,'反三角波');
Xc(9:
32)=0;Xd(9:
32)=0;k=1:
32;figure;
subplot(2,2,1);stem(k-1,Xc);hold;plot(k-1,Xc);
xlabel('n');ylabel('时域特性');text(8,3,'三角波');
subplot(2,2,2);stem(k-1,abs(fft(Xc)));hold;plot(k-1,abs(fft(Xc)));
xlabel('k');ylabel('幅频特性');text(4,10,'三角波');
subplot(2,2,3);stem(k-1,Xd);hold;plot(k-1,Xd);
xlabel('n');ylabel('时域特性');text(2,3,'反三角波');
subplot(2,2,4);stem(k-1,abs(fft(Xd)));hold;plot(k-1,abs(fft(Xd)));
xlabel('k');ylabel('幅频特性');text(4,12,'反三角波');
实验现象分析:
由图知,三角波序列和反三角波序列的时域图像关于横轴成镜像关系,但频域图像完全一样,只是因为幅频图是对x(k)的值取绝对值。
实际频谱相位相反,幅值对应相等。
由实验所得的图形知,N=32点时,两幅度谱都更加密集,更多离散点的幅值显示,“栅栏效应”减小。
在原序列的末端填补零值,变动了DFT的点数,人为的改变了对真实频谱采样的点数和位置,相当于搬动了“尖桩栅栏”的位置,从而使得频谱的峰点和谷点暴露出来。
N=32时,两频谱差别较大,但总体趋势仍然都是中间最小,两侧呈对称。
(4)、一个连续信号含两个频率分量,经采样得
x(n)=sin2π*0.125n+cos2π*(0.125+Δf)nn=0,1……,N-1
已知N=16,Δf分别为1/16和1/64,观察其频谱;当N=128时,Δf不变,其结果有何不同,为什么?
clearall;
N=16;detf=1/16;n=[0:
N-1];
x1(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);
detf=1/64;x2(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);
subplot(2,2,1);stem(n,x1);
xlabel('n');ylabel('时域特性');text(6,1,'N=16,detf=1/16');
subplot(2,2,2);stem(n,abs(fft(x1)));
xlabel('n');ylabel('幅值特性');text(6,4,'N=16,detf=1/16');
subplot(2,2,3);stem(n,x2);
xlabel('n');ylabel('时域特性');text(6,1,'N=16,detf=1/64');
subplot(2,2,4);stem(n,abs(fft(x2)));
xlabel('n');ylabel('幅值特性');text(6,4,'N=16,detf=1/64');
N=128;detf=1/16;n=[0:
N-1];
x3(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);
detf=1/64;
x4(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);
figure;subplot(2,2,1);stem(n,x3);
xlabel('n');ylabel('时域特性');axis([0128-22]);text(6,1.5,'N=128,detf=1/16');
subplot(2,2,2);stem(n,abs(fft(x3)));
xlabel('n');ylabel('幅值特性');axis([0128-1070]);text(40,60,'N=128,detf=1/16');
subplot(2,2,3);stem(n,x3);
xlabel('n');ylabel('时域特性');axis([0128-22]);text(6,1.5,'N=128,detf=1/16');
subplot(2,2,4);stem(n,abs(fft(x4)));
xlabel('n');ylabel('幅值特性');axis([0128-1070]);text(40,60,'N=128,detf=1/16');
实验现象分析:
由图可以看出N=16时,当f由1/16减小为1/64时,频谱图发生较大变化。
此时分辨出1/64对频率的分辨率要求提高,需要通过增大信号序列长度以提高频率的分辨率。
f由1/16减小为1/64时,此时采样是非周期采样,DFT后导致频谱泄漏,栅栏效应使得无法显现出0.125+1/64和0.125两个频率谱峰。
当N增加至128时,频谱更加密集,分辨率提高(1/128),能够分辨出两谱峰。
(5)、用FFT分别实现xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的16点圆周卷积和线性卷积。
clearall;
N=16;
n=0:
N-1;
p=8;q=2;
Xa(n+1)=exp(-(n-p).^2./q);
a=0.1;f=0.0625;
Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n);
Fa=fft(Xa);Fb=fft(Xb);
Fx=Fa.*Fb;
X51=ifft(Fx);
stem(n,X51);
Xa(N+1:
2*N-1)=0;
Xb(N+1:
2*N-1)=0;
Fa=fft(Xa);Fb=fft(Xb);
Fc=Fa.*Fb;
X52=ifft(Fc);
figure;stem(1:
2*N-1,X52);
(6)产生一512点的随机序列xe(n),并用xc(n)和xe(n)作线性卷积,观察卷积前后xe(n)频谱的变化。
要求将xe(n)分成8段,分别采用重叠保留法和重叠相加法。
重叠保留法:
1.在长序列xe(n)前补N1-1=7个零,由于最后输出y(n)要多取8个,为满足L>=N1+N2-1,在xe(n)末尾补8个零,构成长度为527的序列;
2.分段时从xe(n)中每次取72个数,即要做72点的DFT,但每次取数的起始点为每隔65个开始取一次,即相邻两段之间有7个数是相重叠的,也就是将要去掉的部分;
3.将短序列xc(n)作72点的DFT,与每段序列做圆周卷积;
4.去掉所得每段卷积序列的前N1-1=7个点,取其后的65个点。
将各相邻段收尾连结起来,构成长为519的线性卷积序列。
程序代码:
n=1:
512;
xe1=randn(1,512);
n=1:
7;xe(n)=0;
n=8:
519;xe(n)=xe1(n-7);
n=520:
527;xe(n)=0;
n=1:
4;xc(n)=n-1;
n=5:
8;xc(n)=9-n;
n=0:
526
subplot(2,1,1);stem(n,xe);ylabel('xe(k)');title('512点的随机序列');
subplot(2,1,2);n=0:
7;stem(n,xc);
y1=conv(xe1,xc);
figure;subplot(2,1,1);n=0:
518;stem(n,y1);ylabel('y1(n)');title('线性卷积');
hxc=fft(xc,72);
y=zeros(1,527);
fori=0:
7
n=(65*i+1):
65*(i+1)+7;
xi(n-65*i)=xe(n);
hxi=fft(xi,72);
yi=ifft(hxc.*hxi);
k=(65*i+1):
65*(i+1);
y(k)=yi(k-65*i+7);
end
subplot(2,1,2);m=0:
526;stem(m,y);ylabel('y(n)');;title('重叠保留法');
实验结果分析:
实验证明重叠保留法可以用DFT做圆周卷积来表示线性卷积,以减少存储单元和处理时间的延误。
(7)用FFT分别计算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的16点循环相关和线性相关,问一共有多少种结果,他们之间有何异同点。
clearall;
N=16;
n=0:
N-1;
p=8;q=2;
Xa(n+1)=exp(-(n-p).^2./q);
a=0.1;f=0.0625;
Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n);
N=length(Xa);
Fa=fft(Xa,2*N);Fb=fft(Xb,2*N);
Fx=conj(Fa).*Fb;
X71=real(ifft(Fx));
X71=[X71(N+2:
2*N)X71(1:
N)];
n=(-N+1):
(N-1);stem(n,X71);
Xa(N+1:
2*N-1)=0;
Xb(N+1:
2*N-1)=0;
Fa=fft(Xa);Fb=fft(Xb);
Fc=conj(Fa).*Fb;
X72=real(ifft(Fc));
figure;stem(1:
2*N-1,X72);
实验结果分析:
一共可以有16种结果,循环相关是线性相关的一个主值区间,取值的起点不同,其序列即不同,但是一个主值区间包含的信息量是相同的,包含的频率成分都是相同的。
(8)用FFT分别计算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理 吴镇杨 实验