数字信号处理实验快速傅里叶变换及其应用.docx
- 文档编号:7573059
- 上传时间:2023-01-25
- 格式:DOCX
- 页数:27
- 大小:131.71KB
数字信号处理实验快速傅里叶变换及其应用.docx
《数字信号处理实验快速傅里叶变换及其应用.docx》由会员分享,可在线阅读,更多相关《数字信号处理实验快速傅里叶变换及其应用.docx(27页珍藏版)》请在冰豆网上搜索。
数字信号处理实验快速傅里叶变换及其应用
数字信号处理实验报告
第一次实验:
快速傅立叶变换(FFT)及其应用
王宇阳
04011345
一.实验目的:
(1)在理论学习的基础上,通过本实验,加深对FFT的理解,熟悉MATLAB中的有关函数。
(2)应用FFT对典型信号进行频谱分析。
(3)了解应用FFT进行信号频谱分析过程中可能出现的问题,以便在实际中正确应用FFT。
(4)应用FFT实现序列的线性卷积和相关。
二.实验原理:
(1)混叠:
采样序列的频谱是被采样信号频谱的周期延拓,当采样频率不满足奈奎斯特采样定理的时候,就会发生混叠,使得刺痒后的序列信号的频谱不能真实的反映原采样信号的频谱。
(2)泄露:
根据理论分析,一个时间的信号其频带宽度为无限,一个时间无限的信号其频带宽度则为有限。
因此对一个时间有限的信号,应用DFT进行分析,频谱混叠难以避免。
对一个时间无限的信号虽然频带有限,但在实际运算中,时间总是取有限值,在将信号截断的过程中,出现了分散的扩展谱线的现象,称之为频谱泄露或功率泄露。
(3)栅栏效应:
DFT是对单位圆上Z变换的均匀采样,所以它不可能将频谱视为一个连续函数,就在一定意义上看,用DFT来观察频谱就好象通过一个栅栏来观看一个景象一样,只能在离散点上看到真实的频谱,这样就有可能发生一些频谱的峰点和谷点被“尖桩的栅栏”所挡住,不能被我们观察到。
(4)圆周卷积:
把序列X(N)分布在N等份的圆周上,而序列Y(N)经反摺后也分布在另一个具有N等份的同心圆的圆周上。
两圆上对应的数两量两相乘求和,就得到全部卷积序列。
这个卷积过程称做圆周卷积。
(5)互相关函数反映了两个序列X(N)和Y(N)的相似程度,用FFT可以很快的计算互相关函数。
三.实验内容:
实验中用到的函数序列:
(a)Gaussian序列
Xa(n)=exp(-(n-p).^2)/q),0= =0,其他 (b)衰减正弦序列 X(b)=exp(-an)*sin(2pi*fn),0= =0,其他 (c)三角波序列 Xb(n)=n,0= =8-n,4= =0,其他 (d)反三角波序列 Xc(n)=4-n,0= =n-4,4= =0,其他 上机实验内容: 1.观察高斯序列的时域和幅频特性,固定信号xa(n)中参数p=8,改变q的值,使q分别等于2,4,8,观察他们的时域和幅频特性,了解当q取不同值时,对信号序列的时域和幅频特性影响;改变p,使p分别等于8,13,14,观察参数p变化对信号序列的时域和幅频特性影响,注意p等于多少时,会发生明显的泄漏现象,混叠是否也随之出现? 记录实验中观察到的现象,绘出相应的时域序列和幅频特性曲线。 代码: %%高斯序列以p,q为参数生成并输出其时域序列波形以及其dtft; function[timePlot,spectrumPlot]=Gauss(p,q) %%输入相关参数 n=0: 15; N=4000; w=linspace(-pi,pi,N); g=exp(-(n-p).^2/q); spec=dtft(N,g); subplot(2,1,1); stem(n,g);title(strcat('当p=',num2str(p),'q=',num2str(q),'时的高斯时域波形')); subplot(2,1,2); plot(w,spec);title(strcat('当p=',num2str(p),'q=',num2str(q),'时的高斯频域波形')); end %%DTFT变换得到频谱: function[spectrum]=dtft(n,x) %%取4000采样点; w=linspace(-pi,pi,n); spectrum=zeros(1,n); fork=1: n form=1: length(x) wn(k,m)=exp(-1j*w(k)*m); end spectrum(k)=abs(x*wn(k,: )'); end end 实验结果如下: . P=8,q=2 P=8,q=4; P=8,q=8; P=8,q=8; P=13,q=8 P=14,q=8; 分析: 当p固定,q变化时,由时域波形可知,q代表了波形以p为中心向两侧衰减的速度,发现当q越大时,衰减的速度越快,反之,衰减的速度相对较慢,高频分量越小,从频域波形分析可知,三种情况下的形状是相似的,但是三种情况下的频谱在q比较的大的时候在高频上出现了较少的高频分量以及较快的高频分量的衰减。 当q固定,p变化时,由时域波形可知,p代表了波形的移序,在频谱上,显而易见,在p更大的时候产生了更多的高频分量。 右移之后,序列相当于被截断,那么相对应的其频谱泄漏的情况更加的严重,从而出现了高频分量较多的情况。 2、观察衰减正弦序列xb(n)的时域和幅频特性,a=0.1,f=0.0625,检查谱峰出现位置是否正确,注意频谱的形状,绘出幅频特性曲线,改变f,使f分别等于0.4375和0.5625,观察这两种情况下,频谱的形状和谱峰出现的位置,有无混叠和泄漏现象? 说明产生现象的原因。 代码: %%按照alpha和f的值生成相应的正弦序列,画出其时域和频域(dtft)的波形; function[output_args]=sinusoidal(alpha,f) %%生成序列 n=0: 15; N=4000; x=exp(-alpha*n).*sin(2*pi*f*n); %%对离散序列进行dtft变换 spec=dtft(N,x); subplot(2,1,1); stem(n,x);title(strcat('当α=',num2str(alpha),'f=',num2str(f),'时的正弦时域波形')); subplot(2,1,2); w=linspace(-pi,pi,N); plot(w,spec);title(strcat('当α=',num2str(alpha),'f=',num2str(f),'时的正弦频域波形')); end a=0.1,f=0.0625; F=0.4375 F=0.5625; 分析: 在a=0.1固定的条件下,f=0.0625时,谱峰在低频出约等于2*pi*f,但是在f=0.4375和f=0.5626的时候,谱峰主要集中在高频位置。 后两种情况由于不满足nyquist采样定律,所以i出现了混叠。 3.代码; function[]=triangle(N) %UNTITLED5Summaryofthisfunctiongoeshere %Detailedexplanationgoeshere fork=0: 3 x(k+1)=k; end fork=4: 7 x(k+1)=8-k; end spec=fft(x);spec=abs(spec); n=0: N-1; subplot(1,2,1); stem(n,x);title('序列波形'); subplot(1,2,2); stem(n,spec);title(strcat(num2str(N),'点FFT')); End 逆三角波: function[]=reversetriangle(N) %UNTITLED5Summaryofthisfunctiongoeshere %Detailedexplanationgoeshere fork=0: 3 x(k+1)=4-k; end fork=4: 7 x(k+1)=k-4; end fork=8: N-1 x(k+1)=0; end n=0: N-1; spec=fft(x);spec=abs(spec); subplot(2,1,1); stem(n,x);title('序列波形'); subplot(2,1,2); stem(n,spec);title(strcat(num2str(N),'点FFT')); end 正三角波 逆三角波: N=32; 逆三角波 正三角波 分析: 正三角波和逆三角波的时域波形上的走势是不一样的。 在频域上可以发现,N=8时,两者的fft互为相反数,因为两者的和实际上就是一个矩形窗函数,rectwin(8),由于rectwin(8)的fft为80000000,所以除了fft(0)以外,其他都为相反数。 而当N=16时,两者之间的和为一个补零了的rectwin(8), 对于N不同时得到的结果发现,实质上同一个序列补零以后的频谱并没有发生太大的变化,补零后的序列相当了有了更多的采样点,在对应位置处的fft值还是相同的,不同的是,当N变大时,补零长度变大时,暴露出来的更多的采样点值可以减轻fft的栅栏效应。 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不变,其结果有何不同? 代码: function[output_args]=test4(delta_f,n) %UNTITLED7Summaryofthisfunctiongoeshere %Detailedexplanationgoeshere x=zeros(1,n); fork=0: n-1 x(k+1)=sin(2*pi*0.125*k)+cos(2*pi*(0.125+delta_f)*k); end; spec=fft(x); stem(spec); title(strcat('deltaf=',num2str(delta_f),'N=',num2str(n))); end N=16,deltaf=1/16; N=16,deltaf=1/64; N=128,deltaf=1/16 N=128,deltaf=1/64 分析: 由fft频谱图可知,在N=16时,deltaf=1/16时,fft频谱是正确的应该是两个冲激,但是在deltaf=1/64时,此时的信号的带宽减小,由于采样点数是固定的,导致卷积的过程中出现了泄露的情况,因此在第二种情况下出现了谱线扩散的情况。 在N=128时,相当于扩展了加窗的宽度,因而没有出现频谱的泄漏。 5、用FFT分别计算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的16点循环卷积和线形卷积。 代码: n=0: 15; alpha=0.1; f=0.0625; p=8; q=2; k=0: 31; g=exp(-(n-p).^2/q); x=exp(-alpha*n).*sin(2*pi*f*n); fori=17: 32 g(i)=0; x(i)=0; end specg=fft(g); specx=fft(x); convSpec=specg.*specx; linearConv=ifft(convSpec); stem(k,linearConv); 循环卷积: 代码: n=0: 15; alpha=0.1; f=0.0625; p=8; q=2; x=exp(-alpha*n).*sin(2*pi*f*n); g=exp(-(n-p).^2/q); specg=fft(g); specx=fft(x); shiftSpec=specg.*specx; shiftConv=ifft(shiftSpec); stem(n,shiftConv); 6、 代码: clearall; clc; N=512; kn=1: N; a=1: 2*N-1; %%生成512个均值为1,方差为1的随机序列; x=1+1.*randn(512,1); x=x'; fork=1: 4 t(k)=k; end fork=5: 8 t(k)=8-k; end %%直接线性卷积 convout=conv(x,t); specx=fft(t); subplot(3,1,1); stem(convout);title('直接线性卷积'); out=zeros(1,519); out1=zeros(1,519); %%重叠相加法 %%重叠相加法 form=9: 71 t(m)=0; end spec=fft(t); %%重叠相加法 outm=zeros(1,519); fork=1: 8 over(k,: )=[x((64*(k-1)+1): 64*k)zeros(1,7)]; spectrumn(k,: )=fft(over(k,: )).*spec; spectrumn(k,: )=ifft(spectrumn(k,: )); changespectrum(k,: )=[zeros(1,(k-1)*64)spectrumn(k,: )zeros(1,519-64*k-7)]; outm=outm+changespectrum(k,: ); end subplot(3,1,2); stem(outm);title('重叠相加法线性卷积'); %%重叠保留法 xx=[zeros(1,7),x]; fork=1: 8 over1(k,: )=xx(64*(k-1)+1: 64*(k-1)+71); end over1(9,: )=[xx(513: 519)zeros(1,64)]; fork=1: 9 speco1(k,: )=fft(over1(k,: )); spectrumout1(k,: )=speco1(k,: ).*spec; spectrumout1(k,: )=ifft(spectrumout1(k,: )); end s=spectrumout1; out1=[s(1,8: end)s(2,8: end)s(3,8: end)s(4,8: end)s(5,8: end)s(6,8: end)s(7,8: end)s(8,8: end)s(9,8: 14)]; subplot(3,1,3); stem(out1);title('重叠保留法线性卷积'); 7、线性相关: 线性相关即为互相关,在此不再重述。 如题8; 循环相关: clear; p=8;q=2; n=0: 15; xa=exp((-(n-p).^2)/q); a=0.1;f=0.0625; xb=exp(-a*n).*sin(2*pi*f*n); xak=fft(xa); xbk=fft(xb); subplot(2,1,1); yn=ifft(conj(xak).*xbk); stem(n,yn); xlabel('n');ylabel('幅度'); title('循环相关'); subplot(2,1,2); yn=ifft(conj(xbk).*xak); stem(n,yn); xlabel('n');ylabel('幅度'); title('循环相关'); 8.互相关 代码: rxy clear; clc; n=16; alpha=0.1; f=0.0625; p=8; q=2; fori=1: n g(i)=exp(-(i-p)^2/q); x(i)=exp(-alpha*(i-1))*sin(2*pi*f*(i-1)); end fori=n+1: 2*n g(i)=0; x(i)=0; end specx=fft(x,2*n); specg=fft(g,2*n); corr=ifft(conj(specx).*specg); rm=real(corr); rm=[rm(n+2: 2*n)rm(1: n)]; m=(-n+1): (n-1); stem(m,rm); 结果: Ryx clear; clc; n=16; alpha=0.1; f=0.0625; p=8; q=2; fori=1: n g(i)=exp(-(i-p)^2/q); x(i)=exp(-alpha*(i-1))*sin(2*pi*f*(i-1)); end fori=n+1: 2*n g(i)=0; x(i)=0; end specx=fft(x,2*n); specg=fft(g,2*n); corr=ifft(specx.*conj(specg)); rm=real(corr); rm=[rm(n+2: 2*n)rm(1: n)]; m=(-n+1): (n-1); stem(m,rm); 自相关、xa(n) clear; clc; n=16; alpha=0.1; f=0.0625; p=8; q=2; fori=1: n g(i)=exp(-alpha*(i-1))*sin(2*pi*f*(i-1)); x(i)=exp(-alpha*(i-1))*sin(2*pi*f*(i-1)); end fori=n+1: 2*n g(i)=0; x(i)=0; end specx=fft(x,2*n); specg=fft(g,2*n); corr=ifft(conj(specx).*specg); rm=real(corr); rm=[rm(n+2: 2*n)rm(1: n)]; m=(-n+1): (n-1); stem(m,rm); 结果: Xb(n)自相关 代码: clear; clc; n=16; alpha=0.1; f=0.0625; p=8; q=2; fori=1: n g(i)=exp(-(i-p)^2/q); x(i)=exp(-(i-p)^2/q); end fori=n+1: 2*n g(i)=0; x(i)=0; end specx=fft(x,2*n); specg=fft(g,2*n); corr=ifft(conj(specx).*specg); rm=real(corr); rm=[rm(n+2: 2*n)rm(1: n)]; m=(-n+1): (n-1); stem(m,rm); 结果: 三、思考题: 1、实验中的信号序列Xc(n)和Xd(n),在单位圆上的变化频谱|Xc(exp(jw)|和|Xd(exp(jw))|会相同吗? 如果不同,说出哪一个低频分量更多一些,为什么? 答: 由三角波的32点fft可以看出,两者的频谱是不一样的。 逆三角波的频谱的低频分量更多一点。 2、对于一个有限长度序列进行DFT等价于将该序列周期延拓后进行DFT展开,因为DFS也知识取其中一个周期来运算,所以FFT在一定条件下也可以用以分析周期信号序列。 如果实正弦信号sin(2pi*fn)f=0.1用16点FFT来做DFS运算,得到的频谱是信号本身的真实谱吗? 为什么? 答: 不是。 因为原信号的周期是10,因而进行16点FFT时,其16点的周期延拓后的波形已经不再是原来的波形了。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号处理实验 快速傅里叶变换及其应用 数字信号 处理 实验 快速 傅里叶变换 及其 应用