数字信号处理实验快速傅里叶变换及其应用

1、数字信号处理实验报告第一次实验：快速傅立叶变换（FFT）及其应用王宇阳04011345一 实验目的：（1） 在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。（2） 应用FFT对典型信号进行频谱分析。（3） 了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。（4） 应用FFT实现序列的线性卷积和相关。二 实验原理：（1） 混叠：采样序列的频谱是被采样信号频谱的周期延拓，当采样频率不满足奈奎斯特采样定理的时候，就会发生混叠，使得刺痒后的序列信号的频谱不能真实的反映原采样信号的频谱。（2） 泄露：根据理论分析，一个时间的信号其频带宽度为

2、无限，一个时间无限的信号其频带宽度则为有限。因此对一个时间有限的信号，应用DFT进行分析，频谱混叠难以避免。对一个时间无限的信号虽然频带有限，但在实际运算中，时间总是取有限值，在将信号截断的过程中，出现了分散的扩展谱线的现象，称之为频谱泄露或功率泄露。（3） 栅栏效应：DFT是对单位圆上Z变换的均匀采样，所以它不可能将频谱视为一个连续函数，就在一定意义上看，用DFT来观察频谱就好象通过一个栅栏来观看一个景象一样，只能在离散点上看到真实的频谱，这样就有可能发生一些频谱的峰点和谷点被“尖桩的栅栏”所挡住，不能被我们观察到。（4） 圆周卷积：把序列X（N）分布在N等份的圆周上，而序列Y（N）经反摺后

3、也分布在另一个具有N等份的同心圆的圆周上。两圆上对应的数两量两相乘求和，就得到全部卷积序列。这个卷积过程称做圆周卷积。（5） 互相关函数反映了两个序列X（N）和Y（N）的相似程度，用FFT可以很快的计算互相关函数。三 实验内容：实验中用到的函数序列：(a) Gaussian序列Xa(n)=exp(-(n-p).2)/q),0=<n=<15 =0,其他 （b）衰减正弦序列 X(b)=exp(-an)*sin(2pi*fn),0=<n=<15 =0,其他 (c)三角波序列 Xb(n)=n,0=<n=<3 =8-n,4=<n=<7 =0,其他 (d)反

4、三角波序列 Xc(n)=4-n,0=<n=<3 =n-4,4=<n=<7 =0,其他上机实验内容：1 观察高斯序列的时域和幅频特性，固定信号xa(n)中参数p=8，改变q的 值，使q分别等于2，4，8，观察他们的时域和幅频特性，了解当q取不同值时，对信号序列的时域和幅频特性影响；改变p，使p分别等于8，13，14，观察参数p变化对信号序列的时域和幅频特性影响，注意p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。代码：%高斯序列以p，q为参数生成并输出其时域序列波形以及其dtft；function time

5、Plot,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=',num2st

6、r(q),'时的高斯频域波形');end%DTFT变换得到频谱：function spectrum = dtft( n,x )%取4000采样点；w=linspace(-pi,pi,n);spectrum=zeros(1,n);for k=1:n for m=1:length(x) wn(k,m)=exp(-1j*w(k)*m); end spectrum(k)=abs(x*wn(k,:)');end end实验结果如下：.P=8,q=2P=8,q=4;P=8,q=8;P=8,q=8;P=13,q=8P=14,q=8;分析：当p固定，q变化时，由时域波形可知，q代表了波

7、形以p为中心向两侧衰减的速度，发现当q越大时，衰减的速度越快，反之，衰减的速度相对较慢，高频分量越小，从频域波形分析可知，三种情况下的形状是相似的，但是三种情况下的频谱在q比较的大的时候在高频上出现了较少的高频分量以及较快的高频分量的衰减。当q固定，p变化时，由时域波形可知，p代表了波形的移序，在频谱上，显而 易见，在p更大的时候产生了更多的高频分量。右移之后，序列相当于被截断， 那么相对应的其频谱泄漏的情况更加的严重，从而出现了高频分量较多的情况。2、观察衰减正弦序列xb(n)的时域和幅频特性，a=0.1,f=0.0625,检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，

8、使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=

9、',num2str(f),'时的正弦时域波形');subplot(2,1,2);w=linspace(-pi,pi,N);plot(w,spec);title(strcat('当=',num2str(alpha),'f=',num2str(f),'时的正弦频域波形'); enda=0.1,f=0.0625;F=0.4375F=0.5625;分析：在a=0.1固定的条件下，f=0.0625时，谱峰在低频出约等于2*pi*f，但是在f=0.4375和f=0.5626的时候，谱峰主要集中在高频位置。后两种情况由于不满足nyqui

10、st采样定律，所以i出现了混叠。3.代码;function = triangle( N )%UNTITLED5 Summary of this function goes here% Detailed explanation goes herefor k=0:3 x(k+1)=k;endfor k=4:7 x(k+1)=8-k;endspec=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

11、),'点FFT');End逆三角波： function = reversetriangle( N )%UNTITLED5 Summary of this function goes here% Detailed explanation goes herefor k=0:3 x(k+1)=4-k;endfor k=4:7 x(k+1)=k-4;endfor k=8:N-1 x(k+1)=0;endn=0:N-1;spec=fft(x);spec=abs(spec);subplot(2,1,1);stem(n,x);title('序列波形');subplot(2,1

12、,2);stem(n,spec);title(strcat(num2str(N),'点FFT'); end正三角波逆三角波：N=32;逆三角波正三角波分析：正三角波和逆三角波的时域波形上的走势是不一样的。在频域上可以发现，N=8时，两者的fft互为相反数，因为两者的和实际上就是一个矩形窗函数，rectwin(8)，由于rectwin（8）的fft为8 0 0 0 0 0 0 0，所以除了fft（0）以外，其他都为相反数。而当N=16时，两者之间的和为一个补零了的rectwin（8），对于N不同时得到的结果发现，实质上同一个序列补零以后的频谱并没有发生太大的变化，补零后的序列相当

13、了有了更多的采样点，在对应位置处的fft值还是相同的，不同的是，当N变大时，补零长度变大时，暴露出来的更多的采样点值可以减轻fft的栅栏效应。4、一个连续信号含两个频率分量，经采样得x(n)=sin2*0.125n+cos2*(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)%UNTITLED7 Summary of this function goes here% Detailed explanation goes here

14、x=zeros(1,n);for k=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/16N=128,deltaf=1/64分析：由fft频谱图可知，在N=16时，deltaf=1/16时，fft频谱是正确的应该是两个冲激，

15、但是在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);for i=17:32 g(i)=0; x(i)=0;endspecg=ff

16、t(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、代码：clear all;clc;N=512;kn=1:

17、N;a=1:2*N-1;%生成512个均值为1，方差为1的随机序列；x = 1+1.*randn(512,1);x=x'for k=1:4 t(k)=k;endfor k=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);%重叠相加法%重叠相加法for m=9:71 t(m)=0;endspec=fft(t);%重叠相加法outm=zeros(1,519)

18、;for k=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=zero

19、s(1,7),x;for k=1:8 over1(k,:)=xx(64*(k-1)+1:64*(k-1)+71);end over1(9,:)=xx(513:519) zeros(1,64);for k=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:en

20、d) 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'

21、);ylabel('幅度');title('循环相关');subplot(2,1,2);yn=ifft(conj(xbk).*xak);stem(n,yn);xlabel('n');ylabel('幅度');title('循环相关');8. 互相关代码：rxyclear;clc;n=16;alpha=0.1;f=0.0625;p=8;q=2;for i =1:n g(i)=exp(-(i-p)2/q); x(i)=exp(-alpha*(i-1)*sin(2*pi*f*(i-1);endfor i=n+1:2*n

22、g(i)=0; x(i)=0;endspecx=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);结果：Ryxclear;clc;n=16;alpha=0.1;f=0.0625;p=8;q=2;for i =1:n g(i)=exp(-(i-p)2/q); x(i)=exp(-alpha*(i-1)*sin(2*pi*f*(i-1);endfor i=n+1:2*n g(i)=0; x(i)=0;endspecx=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;for i =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);en