仲易+数字信号处理实验报告

北京理工大学2011－2012学年第一学期数字信号处理实验报告1数字信号处理实验报告姓名：仲易班级：05110901学号：20091087北京理工大学2011－2012学年第一学期数字信号处理实验报告2实验1基2-FFT算法实现一、实验目的1.掌握基2-FFT的原理及具体实现方法。2.编程实现基2-FFT算法。3.深刻理解FFT算法的特点。二、实验设备与环境计算机，matlab软件环境。三、实验基础理论FFT是DFT的一种快速算法，能使DFT的计算大大简化，运算时间缩短。FFT利用了WNnk的三个固有特性。即对称性，周期性和可约性，将长序列的DFT分解为短序列的DFT，合并了DFT运算中的某些项，从而减少了DFT的运算量。FFT算法基本上可分为两大类，即按时间抽取法和按频率抽取法。在实现FFT算法时，要重点考虑两个问题，注意数据的读取和存储：(1)输入输出的排序；（2）蝶形运算的实现。按时间抽取算法中输入反序输出顺序，按频率抽取算法中输入顺序输出反序；运算过程中的每一级都有N/2个蝶形运算构成，每一个蝶形运算单元中，两个节点变量运算后得到的结果为下一列相同位置的节点变量，而和其他节点变量无关，可以采用原位运算，节省存储单元。另外，蝶形运算中的复系数WNnk可以存储为能及时查阅的系数表，这样可以借阅运算量，但是需要N/2哥复数存储器。MATLAB中提供了用于计算FFT的函数fft，可将实验中所得到的结果与利用MATLAB中fft函数计算的结果相比较，以此验证结果的正确性。四.实验内容及实验过程1.编程实现序列长度为N=8的按时间抽取的基2-FFT算法。给定一个8点序列，采用编写的程序计算其DFT，并与MATLAB中fft函数计算的结果相比较，以验证结果的正确性。代码如下：x=[12345678];m=3;%X序列长度为2的3次幂N=8;nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1;%求倒序列y=x(nxd);formm=1:m%对x做m次基2分解Nz=2^mm;u=1;WN=exp(-i*2*pi/Nz);forj=1:Nz/2fork=j:Nz:Nkp=k+Nz/2;t=y(kp)*u;%蝶形运算乘项y(kp)=y(k)-t;%加项y(k)=y(k)+t;endu=u*WN;endendy北京理工大学2011－2012学年第一学期数字信号处理实验报告3y1=fft(x)%对比matlab的ffty=Columns1through336.0000-4.0000+9.6569i-4.0000+4.0000iColumns4through6-4.0000+1.6569i-4.0000-4.0000-1.6569iColumns7through8-4.0000-4.0000i-4.0000-9.6569iy1=Columns1through336.0000-4.0000+9.6569i-4.0000+4.0000iColumns4through6-4.0000+1.6569i-4.0000-4.0000-1.6569iColumns7through8-4.0000-4.0000i-4.0000-9.6569i2、编程实现序列长度为N=8的按频率抽取的基2-FFT算法。给定一个8点序列，采用编写的程序计算其DFT，并与MATLAB中fft函数计算的结果相比较，以验证结果的正确性。代码如下：L=2048;M=256;x=ones(1,L);n=0:M-1;h=cos(0.2*pi*n);ticy1=conv(x,h);tocticX=fft(x,L+M-1);北京理工大学2011－2012学年第一学期数字信号处理实验报告4H=fft(h,L+M-1);Y=X.*H;y2=ifft(Y);toc运行结果如下：（分别为线形卷积和快速卷积）Elapsedtimeis0.003352seconds.Elapsedtimeis0.007167seconds.L=4096;M=256;x=ones(1,L);n=0:M-1;h=cos(0.2*pi*n);ticy1=conv(x,h);tocticX=fft(x,L+M-1);H=fft(h,L+M-1);Y=X.*H;y2=ifft(Y);toc运行结果如下：（分别为线形卷积和快速卷积）Elapsedtimeis0.003666seconds.Elapsedtimeis0.030134seconds.3.将上述FFT程序推广大序列长度为N=2v的情况，要求利用原位运算。代码如下：L=2048;M=256;n=0:M-1;h=cos(0.2*pi*n);l=128;x=ones(1,l);q=0;ticH=fft(h,l+M-1);X=fft(x,l+M-1);Y=X.*H;y=ifft(Y);yp=y;fori=2:16q=q+l;yq=[zeros(1,q),y];北京理工大学2011－2012学年第一学期数字信号处理实验报告5yp=[yp,zeros(1,l)];yp=yq+yp;endtoc运行结果如下：Elapsedtimeis0.018875seconds.4.L=2048;M=256;n=0:M-1;h=cos(0.2*pi*n);x=ones(1,L);Lx=2048;M=256;M1=M-1;N=512;L=N-M1;tich=[h,zeros(1,N-M)];x=[zeros(1,M1),x,zeros(1,N-1)];a=floor((Lx+M1-1)/(L))+1;Y=zeros(1,N);fork=0:a-1xk=x(k*L+1:k*L+N);b=fft(xk,N);C=fft(h,N);Z=b.*C;Y(k+1,:)=ifft(Z,N);endY=Y(:,M:N)';Y=(Y(:))'toc运行结果如下：Elapsedtimeis0.035329seconds.北京理工大学2011－2012学年第一学期数字信号处理实验报告6实验2利用DFT分析信号频谱一、实验目的1、加深对DFT原理的理解。2、应用DFT分析信号的频谱。3、深刻理解利用DFT分析信号频谱的原理，分析实现过程中出现的现象及解决方法。二、实验设备与环境计算机、MATLAB软件环境三、实验基础理论1、DFT与DTFT的关系序列()xn的N点DFT()Xk实际上就是()xn序列的DTFT在N个等间隔频率点2/(01)kkNkN上的采样样本。2、利用DFT求DTFT方法一：利用差值运算12()()()NjkkxeXkN其中()x为内插函数12sin(/2)()sin(/2)NjNeN方法二：在实际MATLAB计算中，采用增加数据的长度N，使得到的DFT谱线更加精细，用其包络近似计算DTFT。如果数据量不足，可采用补零来增加数据长度。3、利用DFT分析连续时间信号的频谱1）确定时域采样间隔T，得到离散序列()xn2）确定截取长度M，得到M点离散序列()()()Mxnxnwn，这里()wn为窗函数3）确定频域采样点数N，要求NM4）利用FFT计算离散序列的N点DFT，得到()MXk5）由()MXk计算()aXj采样点的近似值。2120()()()MjknNaaMknNTXjTxnTeTXk4、可能用到的MATLAB函数与代码实验中DFT运算采用MATLAB中提供的函数fft来实现。DTFT可以利用MATLAB矩阵运算的方法进行计算。北京理工大学2011－2012学年第一学期数字信号处理实验报告712112()[][[],[],,[]]NNjnjnnjjnNnnjneeXexnexnxnxne四、实验内容1、（1）n=0:3;x=[2-111];w=-pi:0.01*pi:pi;X=x*exp(-j*n'*w);subplot(211)plot(w,abs(X));axistightsubplot(212)plot(w,angle(X)/pi);axistight(2)在原基础上增加代码：n=0:3;x=[2-111];N=4;北京理工大学2011－2012学年第一学期数字信号处理实验报告8k=0:3;Y=fft(x);w=0:0.01*pi:2*pi;X=x*exp(-j*n'*w);m=0:2*pi/4:3/2*pisubplot(211);plot(w,abs(X));xlabel('\Omega/pi');title('DTFT-Magnitude');axistightholdonstem(m,abs(Y));subplot(212);plot(w,angle(X)/pi);xlabel('\Omega/pi');title('DTFT-Phase');holdonstem(m,angle(Y))（3）在原本基础上设置代码subplot(211);Z=fft(x,64);%n=64DFTN=0:63;stem(abs(Z));subplot(212);stem(angle(Z))北京理工大学2011－2012学年第一学期数字信号处理实验报告92、考察序列()cos(0.48)cos(0.52)xnnn（1）010n，用DFT估计()xn的频谱；将()xn补零加长到长度为100点序列用DFT估计()xn的频谱。要求画出相应波形。程序代码：程序输出结果：北京理工大学2011－2012学年第一学期数字信号处理实验报告10补零后：程序代码：程序输出结果：北京理工大学2011－2012学年第一学期数字信号处理实验报告11（2）0100n时，用DFT估计()xn的频谱，并画出波形。程序代码：程序输出结果：北京理工大学2011－2012学年第一学期数字信号处理实验报告12（3）根据实验结果，分析怎样提高频谱分辨率。答：减小截取长度或减小采样间隔的方法可以提高频谱分辨率。3、已知信号123()0.15sin(2)sin(2)0.1sin(2),xtftftft其中1f=1Hz，2f=2Hz，3f=3Hz。从()xt的表达式可以看出，它包含三个频率的正弦波，但是，从其时域波形来看，似乎是一个正弦信号，利用DFT做频谱分析，确定适合的参数，使得到的频谱的频率分辨率符合需要。程序代码：程序输出结果：北京理工大学2011－2012学年第一学期数字信号处理实验报告13利用DFT近似分析连续时间信号0.1()()txteut的频谱。分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定适合的参数。程序代码：程序输出结果：北京理工大学2011－2012学年第一学期数字信号处理实验报告14通过改变步长可以达到改变采样间隔的结果，输出不同的频谱。北京理工大学2011－2012学年第一学期数字信号处理实验报告15实验5、脉冲响应不变法设计IIR数字滤波器一、实验目的1、掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法2、加深理解数字滤波器和模拟滤波器之间的技术指标转化。3、掌握脉冲相应不变法设计IIR数字滤波器的优缺点及适用范围。二、实验设备与环境计算机、MATLAB软件环境三、实验基础理论1、基本原理：从时域响应出发，使数字滤波器的单位脉冲响应h（n）模仿模拟滤波器的单位冲激响应。2、设计步骤：a.确定数字滤波器性能指标。b.将数字滤波器频率指标转换成相应的模拟滤波器的频率指标c.设计模拟滤波器实现方法：脉冲响应不变法[r,p,k]=residue(b,a)[b,a]=residue(r,p,k)d.模拟滤波器数字化实现方法：[bz,az]=impinvar(b,a,fs)四、实验内容1、设采样频率fs=4khz,采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其3db截止频率为fc=1khz.首先设计巴特沃斯滤波器：N=3;Wc=0.25*pi;[b,a]=butter