数字信号处理实验报告.doc

完美.格式.编辑 数字信号处理实验讲义 专业.资料.整理 完美.格式.编辑 前 言2实验一 MATLAB简介3实验二 用FFT实现信号的谱分析5实验三 IIR数字巴特沃思滤波器的设计8实验四 FIR数字滤波器的设计9 专业.资料.整理 实验一 MATLAB简介实验目的1熟悉MATLAB软件的使用方法；2MATLAB的绘图功能；3用MATLAB语句实现信号的描述及变换。实验原理1在MATLAB下编辑和运行程序在MATLAB中，对于简单问题可以在命令窗（command windows）直接输入命令，得到结果；对于比较复杂的问题则可以将多个命令放在一个脚本文件中，这个脚本文件是以m为扩展名的，所以称之为M文件。用M文件进行程序的编辑和运行步骤如下：（1）打开MATLAB，进入其基本界面；（2）在菜单栏的File项中选择新建一个M文件；（3）在M文件编辑窗口编写程序；（4）完成之后，可以在编辑窗口利用Debug工具调试运行程序，在命令窗口查看输出结果；也可以将此文件保存在某个目录中，在MATLAB的基本窗口中的File项中选择Run The Script，然后选择你所要运行的脚本文件及其路径，即可得出结果；也可以将此文件保存在当前目录中，在MATLAB命令窗口，“”提示符后直接输入文件名。2MATLAB的绘图功能plot(x,y) 基本绘图函数，绘制 x和y之间的坐标图。 figure(n ) 开设一个图形窗口n subplot(m,n,N) 分割图形窗口的MATLAB函数，用于在一个窗口中显示多个图形，将图形窗口分为m行n列，在第N个窗口内绘制图形。axis(a0,b0,a1,b1 ) 调整坐标轴状态 title( ) 给图形加题注xlabel ( ) 给x轴加标注ylabel ( ) 给y轴加标注grid 给图形加网格线3信号描述及变换信号描述及变换包括连续时间信号和离散时间信号内容，详细内容请见课本第1章、第2章。实验内容1试用MATLAB绘制出下列信号的波形：(Signal 1.6)（1） ；（2）（3） ；（4） ；（5） 【程序代码】clear all;close all;clc;syms t;x1=exp(-1.5*t)x2=3*sin(0.5*pi*t)x3=0.5+0.5*sym(sign(t)x4=sym(heaviside(t)+sym(heaviside(t-1)-sym(2*heaviside(t-2)x5=.5*t*(sym(heaviside(t)-sym(heaviside(t-4)subplot(2,3,1);ezplot(x1);axis(-6 3 -500 7000);title(x1(t)=exp(-1.5t);grid onsubplot(2,3,2);ezplot(x2);title(x2(t)=3sin(0.5t);grid onsubplot(2,3,3);fplot(sign(t)/2 + 1/2,-10 10,1e-8);ezplot(x3,-10 10);axis(-10 10 -.2 1.2);xlabel(t);title(x3(t)=0.5+0.5sgn(t);grid onsubplot(2,3,4);ezplot(x4,-1 3);title(x4(t)=u(t)+u(t-1)-2u(t-2);grid onsubplot(2,3,5);ezplot(x5,-2 6);title(x5(t)=0.5tu(t)-u(t-4);grid onsubplot(2,3,6);axis off2 已知连续时间信号(Signal 1.7)，试用MATLAB绘制出下列信号的波形：（1） ；（2） ；（3） ；（4） ；（5） 。【程序代码】clear all;close all;clc;figure(2)syms t;x1=(4-t)*(sym(heaviside(t)-sym(heaviside(t-4)x2=exp(-2*t)*sym(heaviside(t)x3=sin(2*pi*t)x4=subs(x1,t,t/2)x5=subs(x4,t,t-2)x6=subs(x2,t,-t)x7=x2+x6x8=x7*x3subplot(2,3,1)text(0,0.9,x1(t)=(4-t)u(t)-u(t-4) );text(0,0.7,x2(t)=exp(-2t)u(t);text(0,0.5,x3(t)=sin(2t);axis off;box off;subplot(2,3,2);ezplot(x5,-3,10);gridtitle(x4(t)=x1(t/2);subplot(2,3,3);ezplot(x5,-1,12);gridtitle(x5(t)=x4(t-2);subplot(2,3,4);ezplot(x6,-10,5);axis(-5 5 -0.1 1.1);gridtitle( x6(t)=x2(-t);subplot(2,3,5);ezplot(x7,-5 5);axis(-5 5 -0.1 1.1);grid;title(x7(t)=x2(t)+x6(t);subplot(2,3,6);tv8=-2.5:0.05:2.5;xv8=subs(x8,tv8);plot(tv8,xv8)axis(-2.5 2.5 -1 1);gridxlabel(t);title(x8(t)=x7(t)x3(t);clear tv8 xv83列出单位冲激信号、单位阶跃信号、正弦信号的MATLAB表达式，并绘出信号波形。【程序代码】clear all;close all;clcsyms t;x1=sym(dirac(t);x2=sym(Heaviside(t);x3=sin(t);tn=-6.3:0.1:6.3;xn1=subs(x1,t,tn);xn2=subs(x2,t,tn);xn3=subs(x3,t,tn);plot(tn,xn1,k,tn,xn2,r,tn,xn3,m);gridxlabel(t);ylabel(x(t);legend(dirac,Heaviside,sin)hold onplot(tn,xn1,k.,tn,xn2,r.,tn,xn3,m.)实验二 用FFT实现信号的谱分析实验目的1了解FFT在信号谱分析中的作用；2了解谱分析的一般步骤和方法。实验原理关于信号谱分析的步骤和方法参见教材第3章相关内容。为了解信号的特点，了解信号频谱分布情况，应该对信号进行谱分析，计算出信号的幅度谱、相位谱和功率谱。信号的谱分析可以用FFT实现，讨论如下：1 谱分析中的参数选择；A 若已知信号的最高频率，为防止混叠，选定采样频率： (1)B 根据实际需要，选定频率分辨，一但选定后，即可确定FFT所需的点数N (2)我们希望越小越好，但越小，N越大，计算量、存储量也随之增大。一般取N为2的整次幂，以便用FFT计算，若已给定N，可用补零方法便N为2的整次幂。C 和N确定后，即可确定所需相应模拟信号的长度 (3)分辨率反比于T，而不是N，在给定的T的情况下，靠减小来增加N是不能提高分辨率的，因为为常数 2谱分析步骤； A 数据准备 (4) B 使用FFT计算信号的频谱 (5) (6) C 由频谱计算幅度谱、相位谱和功率谱 (7) (8) (9)3实验中用到的一些基本函数简介y=fft(x,n) ； 计算n点的FFT。abs(x) ； 取绝对值。angle(z) ； 取相角。Pxx, f= periodogram (xn, nfft, fs, window) ；%周期图谱估计Pxx, f=pwelch (xn, nfft, fs, window, noverlap)；%平均周期图法Pxx=psd (xn) ；功率谱密度 实验内容1. 已知序列x(n)=2sin(0.48n)+cos(0.52n) 0n100,试绘制x(n)及它的频谱图。若x(n)=sin(0.56n)+2cos(0.25n),结果又如何？【程序代码】clear all;close all;clcN=100;n=0:N-1;xn=2*sin(0.48*pi*n)+cos(0.52*pi*n);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1)plot(n,xn)xlabel(n);ylabel(x(n);title(x(n)N=100);subplot(1,2,2)k=0:length(magXK)-1;stem(k,magXK,.);xlabel(k);ylabel(|X(K)|);title(X(K)N=100);2. 对下面信号进行频谱分析，求幅度谱和相位谱。(1)，(2)，【程序代码】clear all;close all;clcfs1=5000;dt1=1/fs1;N1=0.004/dt1;n1=0:N1-1;xn1=0.8.(n1*dt1);Xk1=fft(xn1,N1);mag1=2*abs(Xk1)/N1;pha1=angle(Xk1);f1=n1*fs1/N1;subplot(221)stem(f1,mag1,fill);title(Magnitude 1)xlabel(f);ylabel(|Xk(1)|);subplot(222)stem(f1,pha1,fill);title(Phase 1)xlabel(f);ylabel(|Phk(1)|);fs2=128;dt2=1/fs2;N2=18;n2=0:N2-1;xn2=sin(n2*dt2)./(n2*dt2+eps);Xk2=fft(xn2,N2);mag2=2*abs(Xk2)/N2;pha2=angle(Xk2);f2=n2*fs2/N2;subplot(223)stem(f2,mag2,fill);title(Magnitude 2)xlabel(f);ylabel(|Xk(2)|);subplot(224)stem(f2,pha2,fill);title(Phase 2)xlabel(f);ylabel(|Phk(2)|)3. 给定信号，现在对采样，采样点数，采样频率=50Hz，设采样序列为，编写程序计算的频谱，并绘图；改变采样频率，得到序列，计算的频谱，并绘图；增大采样点数，得到序列，计算的频谱，并绘图；采样点数N=64，采样频率=300Hz，在采样点后补零得到新序列，计算的频谱，并绘图。(Signal 3.18)【程序代码】clear all;close all; clcN=16;n=0:N-1;f1=15;f2=18;fs=50;dt=1/fs;xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;figure(1);subplot(121)plot(f,mag);gridtitle(Magnitude);xlabel(f);ylabel(|Xk|);subplot(122)plot(f,pha);gridtitle(Phase);xlabel(f);ylabel(?Xk?);N=16;n=0:N-1;f1=15;f2=18;fs=36; dt=1/fs;xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;figure(2);subplot(221)plot(f,mag,r);gridtitle(Magnitude(36Hz Sampling requency);xlabel(f);ylabel(|Xk|);subplot(222)plot(f,pha,r);gridtitle(Phase(36Hz Sampling requency);xlabel(f);ylabel(?Xk?);fs=100; dt=1/fs;xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;subplot(223)plot(f,mag,r);gridtitle(Magnitude(100Hz Sampling Frequency);xlabel(f);ylabel(|Xk|);subplot(224)plot(f,pha,r);gridtitle(Phase(100Hz Sampling Frequency);xlabel(f);ylabel(?Xk?);N=64;n=0:N-1;f1=15;f2=18;fs=50;dt=1/fs;xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;figure(3);subplot(221)plot(f,mag,m);gridtitle(Magnitude(64 Sampling Points);xlabel(f);ylabel(|Xk|);subplot(222)plot(f,pha,m);gridtitle(Phase(64 Sampling Points);xlabel(f);ylabel(?Xk?);N=512; n=0:N-1;dt=1/fs;xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;subplot(223)plot(f,mag,m);gridtitle(Magnitude(512 Sampling Points);xlabel(f);ylabel(|Xk|);subplot(224)plot(f,pha,m);gridtitle(Phase(512 Sampling Points);xlabel(f);ylabel(?Xk?);N=64;n=0:N-1;f1=15;f2=18;fs=300;dt=1/fs;xn0=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);xn=xn0,zeros(1,64);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;figure(4);subplot(221)plot(f,mag,g);gridtitle(Magnitude(filled with 64 zeros);xlabel(f);ylabel(|Xk|);subplot(222)plot(f,pha,g);gridtitle(Phase (filled with 64 zeros);xlabel(f);ylabel(?Xk?);xn=xn0,zeros(1,448); Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;subplot(223)plot(f,mag,g);gridtitle(Magnitude(filled with 448 zeros);xlabel(f);ylabel(|Xk|);4. 试求下列差分方程所描述的输出序列的功率谱并作图。(a)(b)(c)式中，是方差为(例如，=1/12)的白噪声。【程序代码】clear all;close all;clcnfft=512;window=hanning(256);noverlap=128;w=randn(1,1000);for i=1:1000if i=2x(i)=randn(1,1);elsex(i)=-0.81*x(i-2)+w(i)-w(i-1);endendPxx f=psd(x,nfft,2,window,noverlap,mean);figure(1)plot(f,10*log10(Pxx);gridxlabel(f);ylabel(Power Spectrum,dB);title(a): x(n)=-0.81x(n-2)+w(n)-w(n-1)w=randn(1,1000);for i=1:1000if i=2x(i)=randn(1,1);elsex(i)=w(i)-w(i-2);endendPxx f=psd(x,nfft,2,window,noverlap,mean);figure(2)plot(f,10*log10(Pxx),m);gridxlabel(f);ylabel(Power Spectrum,dB);title(b): x(n)=w(n)-w(n-2)w=randn(1,1000);for i=1:1000if i=2x(i)=randn(1,1);elsex(i)=-0.81*x(i-2)+w(2);endendPxx f=psd(x,nfft,2,window,noverlap);figure(3)plot(f,10*log10(Pxx),r);gridxlabel(f);ylabel(Power Spectrum,dB);title(c): x(n)=-0.81x(n-2)+w(n)5. 一序列是由两个频率相距为的模拟信号采样得来的，即 n=0,1,15已知序列长度N=16，试采用周期图法，应用DFT分别计算当=0.06及=0.01时的功率谱估计，并通过作图说明从功率谱估计的分布是否能分辨出这两个正弦信号的真实频谱？若N=64又有什么变化？【程序代码】clear all;close all;clcnfft=512;window=hanning(256);noverlap=128;df=0.06 0.01;N=16;n=0:N-1;x1=sin(2*pi*0.135*n)+cos(2*pi*(0.135+df(1)*n);x2=sin(2*pi*0.135*n)+cos(2*pi*(0.135+df(2)*n);Pxx1 f1=psd(x1,nfft,1,window,noverlap,mean);Pxx2 f2=psd(x2,nfft,1,window,noverlap,mean);figure(1)plot(f1,Pxx1);gridxlabel(f);ylabel(Power Spectrum);title(sin2?(0.135)n+cos2?(0.135+0.06)n N=16)figure(2)plot(f2,Pxx2);gridxlabel(f);ylabel(Power Spectrum );title(sin2?(0.135)n+cos2?(0.135+0.01)n N=16)N=64;n=0:N-1;x1=sin(2*pi*0.135*n)+cos(2*pi*(0.135+df(1)*n);x2=sin(2*pi*0.135*n)+cos(2*pi*(0.135+df(2)*n);Pxx1 f1=psd(x1,nfft,1,window,noverlap,mean);Pxx2 f2=psd(x2,nfft,1,window,noverlap,mean);figure(3)plot(f1,Pxx1,m);gridxlabel(f);ylabel(Power Spectrum);title(sin2?(0.135)n+cos2?(0.135+0.06)n N=64)figure(4);plot(f2,Pxx2,m);gridxlabel(f);ylabel(Power Spectrum);title(sin2?(0.135)n+cos2?(0.135+0.01)n N=64)6. 用MATLAB产生256点白噪声序列，应用Welch法估计其功率谱，每段长64点，重叠32点，输出平均后的功率谱图以及对256点一次求周期图的功率谱图。【程序代码】N=256;x=rand(1,N);window=hanning(64);noverlap=32;nfft=128;Pxx1 f1=periodogram(x);figure(1)plot(f1/pi,10*log10(Pxx1);gridxlabel(f);ylabel(Power Spectrum,dB);title(256)Pxx2 f2=psd(x,nfft,2,window,noverlap,line);figure(2)plot(f2,10*log10(Pxx2),m);gridxlabel(f);ylabel(Power Spectrum,dB);title(Welch)7. 离散信号序列，长度N=512，令取值分别为16Hz和15Hz，计算序列的功率谱，比较谱图的差别。试采用不同的窗函数，再比较谱图的变化。【程序代码】clear all;close all;clcN=512;n=0:N-1;nfft=512;noverlap=128;window=hanning(256);fs=32;f=16;x=2*sin(2*pi*f*n/fs);Pxx f=psd(x,nfft,fs,window,noverlap,mean);subplot(321)plot(f,Pxx);gridxlabel(f);ylabel(Power Spectrum);title(hanning window f=16)f=15;x=2*sin(2*pi*f*n/fs);Pxx f=psd(x,nfft,fs,window,noverlap,mean);subplot(322)plot(f,Pxx);gridxlabel(f);ylabel(Power Spectrum);title(hanning window f=15)window=hanning(256);f=16;x=2*sin(2*pi*f*n/fs);Pxx f=psd(x,nfft,fs,window,noverlap,mean);subplot(323)plot(f,Pxx);gridxlabel(f);ylabel(Power Spectrum);title( boxcar window f=16)f=15;x=2*sin(2*pi*f*n/fs);Pxx f=psd(x,nfft,fs,window,noverlap,mean);subplot(324)plot(f,Pxx);gridxlabel(f);ylabel(Power Spectrum);title( boxcar window f=15)window=kaiser(256);f=16;x=2*sin(2*pi*f*n/fs);Pxx f=psd(x,nfft,fs,window,noverlap,mean);subplot(325)plot(f,Pxx);gridxlabel(f);ylabel(Power Spectrum);title( kaiser window f=16)f=15;x=2*sin(2*pi*f*n/fs);Pxx f=psd(x,nfft,fs,window,noverlap,mean);subplot(326)plot(f,Pxx);gridxlabel(f);ylabel(Power Spectrum);title( kaiser window f=15)8. 已知一个被白噪声污染的信号，其中，=25Hz，=75Hz，=150Hz。应用Welch法进行功率谱估计并绘制功率谱图。【程序代码】clear;close;clcN=512;n=0:N-1;f1=25;f2=75;f3=150;Fs=500;Ts=1/Fs;noverlap=80;nfft=256;rn=randn(1,N);window=hanning(128);xn=2*sin(2*pi*f1*n*Ts)+0.5*sin(2*pi*f2*n*Ts)+0.5*sin(2*pi*f3*n*Ts)+rn;Pxx f=psd(xn,nfft,Fs,window,noverlap,mean);plot(f,10*log(Pxx);gridxlabel(f/Hz);ylabel(Power Spectrum,dB);title(Welch f1=25Hz,f2=75Hz,f3=150Hz)figureplot(f,Pxx,m);gridxlabel(f/Hz);ylabel(Power Spectrum);title(Welch f1=25Hz,f2=75Hz,f3=150Hz)实验三 IIR数字巴特沃思滤波器的设计实验目的1 掌握用模拟滤波器原型设计IIR滤波器的基本方法；2 掌握数字巴特沃思滤波器的设计方法与步骤；3 理解系统频率响应的概念，学习编写计算系统频率响应的方法。实验原理1数字巴特沃思滤波器设计的详细内容见教材第4章，现将设计步骤归纳如下：A 根据给定的频带指标，由双线性变换的频率关系，确定相应的模拟滤波器原型频带指标； B 利用原型低通滤波器，选择合适的参数，设计出符合指标的模拟低通滤波器； C 利用双线性变换，将所获得的模拟滤波器的s域表示转换为相应数字滤波器的z域表示，即它的系统函数，再利用IIR滤波器设计方案具体实现该滤波器。 2. 实验中用到的一些基本函数见教材第4.5.1节（与IIR数字滤波器设计相关MATLAB函数）。实验内容1 希望设计一个巴特沃斯低通数字滤波器，其3dB带宽为0.2，阻带边缘频率为0.5，阻带衰减大于30dB。给定采样间隔。用双线性变换法设计该低通数字滤波器。给出它的及对数幅频响应。【程序代码】clear all;close all;clcTs=10e-6;fs=1/Ts;wp=0.2*pi;ws=0.5*pi;rp=3;rs=30;n wn=buttord(wp/pi,ws/pi,rp,rs);numz denz=butter(n,wp/pi);Gz=tf(numz,denz,Ts)h w=freqz(numz,denz,128,fs);plot(w,abs(h),r);xlabel(f);ylabel(H);title(LowPass Filter);grid on2 给定待设计的数字高通和带通滤波器的技术指标如下：(1) HP：，。(2) BP：，。试用双线性变换分别设计满足上述要求的巴特沃斯滤波器，给出其系统函数、对数幅频及相频曲线。【程序代码】clear all;close all;clcFs=1000;Ts=1/Fs;fp=400;fs=300;wp=2*pi*fp*Ts;ws=2*pi*fs*Ts;rp=3;rs=35;n wn=buttord(wp/pi,ws/pi,rp,rs);numz denz=butter(n,wp/pi,high);Gz_HP=tf(numz,denz,Ts)h w=freqz(numz,denz,128,Fs);figure(1)subplot(2,2,1 2)plot(w,20*log10(abs(h),r);gridxlabel(f);ylabel(20log(H);box offtitle(High Pass Filter);subplot(2,2,3)plot(w,abs(h),r);xlabel(f);ylabel(H);box offsubplot(224)plot(w,angle(h),m);gridxlabel(f);ylabel(phase);grid on;box offclear all;Fs=2000;Ts=1/Fs;fs1=200;f1=300;f2=400;fsh=500;wp=2*pi*fs1*Ts;ws=2*pi*fsh*Ts;w1=2*pi*f1*Ts;w2=2*pi*f2*Ts;rp=3;rs=40;n wn=buttord(wp/pi,ws/pi,rp,rs);numz denz=butter(n,w1/pi w2/pi);Gz_BP=tf(numz,denz,Ts)h w=freqz(numz,denz,128,Fs);figure(2)subplot(2,2,1 2);plot(w,20*log10(abs(h),r);xlabel(f);ylabel(20log(H);grid on;box offtitle(Band Pass Filter);subplot(2,2,3);plot(w,abs(h),r);gridxlabel(f);ylabel(H);grid on;box offsubplot(2,2,4);plot(w,angle(h),m);xlabel(f);ylabel(phase);grid on;box offclear all实验四 FIR数字滤波器的设计实验目的1 掌握FIR数字滤波器的设计方法与步骤；2 理解系统频率响应的概念，学习编写计算系统频率响应的方法。实验原理1. FIR数字滤波器设计的详细内容见教材第4章。(1) 窗口法：窗口法设计FIR数字滤波器的步骤： A 给出希望的滤波器频率响应函数； B 根据允许的过渡带宽度及阻带衰减确定所采用的窗函数和N值； C 做的逆傅里叶变换得hd(n)； D对hd(n)加窗处理得到有限长序列 E对h(n)做傅里叶变换得到频率响应，用作为的逼近，并用给定的技术指标来检验。(2) 切比雪夫一致逼近计算机辅助设计方法：用Remez算法实现FIR滤波器的等波纹逼近。详细内容教材第4.4.3节。 2. 实验中用到的一些基本函数见教材第第4.5.2节（与FIR数字滤波器设计相关MATLAB函数）。实验内容1. 请选择合适的窗函数及N来设计一个线性相位低通滤波器 要求其最小阻带衰减为45dB，过渡带宽为，(1) 已知，求出并画出曲线。(2) 保留原有轨迹，画出用满足所给条件的其他几种窗函数设计出的曲线。【程序代码】clear all;close all;clcrs=45;w0=8*pi/51;w0_ham=3.3; N=2*pi*w0_ham/w0;M=ceil(N-1)-1;wc=0.5*pi;hn=fir1(M,wc/pi,hamming(M+1)s