matlab实验指导.doc

数字信号处理MATLAB实验指导实验一：离散时间信号和离散时间系统 一、 实验目的：1、 以MATLAB为工具学习离散时间信号的图形表示；2、 以课本提供的例程，练习、理解典型离散信号及其基本运算；3、 通过MATLAB仿真简单的离散时间系统，研究其时域特性；4、 加深对离散系统的差分方程、冲激响应和卷积分析方法的理解。二、 实验内容：1、 典型序列的产生与显示；2、 序列的简单运算；3、 复合和复杂信号的产生与显示；4、 离散时间系统的仿真：线性和非线性系统、时变和非时变系统的仿真；5、 线性时不变离散时间系统：系统冲激响应、卷积运算、系统的级联、系统的稳定性；三、 实验例程： 1、 参照课本例程产生下列序列，并用plot、stem好人subplot命令绘出其图形：单位取样序列；单位阶跃序列；矩形序列RN(n)，斜变序列。所需输入的数据是产生序列的长度L和抽样频率FT。% Program P1_1% Generation of a Unit Sample Sequence clf;% Generate a vector from -10 to 20n = -10:20;% Generate the unit sample sequenceu = zeros(1,10) 1 zeros(1,20);% Plot the unit sample sequencestem(n,u);xlabel(Time index n);ylabel(Amplitude);title(Unit Sample Sequence);axis(-10 20 0 1.2);2、 编写MATLAB实指数序列程序，% Program P1_3% Generation of a real exponential sequenceclf;n = 0:35; a = 1.2; K = 0.2;x = K*a.n;stem(n,x);xlabel(Time index n);ylabel(Amplitude);3、 编写产生swept frequency sinusoidal 序列的程序。% Program P1_7% Generation of a swept frequency sinusoidal sequencen = 0:100;a = pi/2/100;b = 0;arg = a*n.*n + b*n;x = cos(arg);clf;stem(n, x);axis(0,100,-1.5,1.5);title(Swept-Frequency Sinusoidal Signal);xlabel(Time index n);ylabel(Amplitude);grid; axis; 4、 产生正弦振幅调制序列% Generation of amplitude modulated sequenceclf;n = 0:100; m = 0.4;fH = 0.1; fL = 0.01; xH = sin(2*pi*fH*n);xL = sin(2*pi*fL*n);y = (1+m*xL).*xH;stem(n,y);grid;xlabel(Time index n);ylabel(Amplitude);5、 用滑动平均滤波器平滑带噪信号，讨论滤波器长度对平滑效果、输出平滑后信号与输入带噪信号之间延时的影响。% Program P1_5% Signal Smoothing by Averagingclf;R = 51;d = 0.8*(rand(R,1) - 0.5); % Generate random noisem = 0:R-1;s = 2*m.*(0.9.m); % Generate uncorrupted signalx = s + d; % Generate noise corrupted signalsubplot(2,1,1);plot(m,d,r-,m,s,g-,m,x,b-.);xlabel(Time index n);ylabel(Amplitude);legend(dn ,sn ,xn );x1 = 0 0 x;x2 = 0 x 0;x3 = x 0 0;y = (x1 + x2 + x3)/3;subplot(2,1,2);plot(m,y(2:R+1),r-,m,s,g-);legend( yn ,sn );xlabel(Time index n);ylabel(Amplitude);6、编写输入序列、计算输出序列、差分输出并画出输出序列。% Program P2_4% Generate the input sequencesclf;n = 0:40; D = 10;a = 3.0;b = -2;x = a*cos(2*pi*0.1*n) + b*cos(2*pi*0.4*n);xd = zeros(1,D) x;num = 2.2403 2.4908 2.2403;den = 1 -0.4 0.75;ic = 0 0; % Set initial conditions% Compute the output yny = filter(num,den,x,ic);% Compute the output ydnyd = filter(num,den,xd,ic);% Compute the difference output dnd = y - yd(1+D:41+D);% Plot the outputssubplot(3,1,1)stem(n,y);ylabel(Amplitude); title(Output yn); grid;subplot(3,1,2)stem(n,yd(1:41);ylabel(Amplitude);title(Output due to Delayed Input xn ?, num2str(D),); grid;subplot(3,1,3)stem(n,d);xlabel(Time index n); ylabel(Amplitude);title(Difference Signal); grid;7、 编写输入序列和系统单位脉冲响应的卷积程序并画出图形。% Program P2_7clf;h = 3 2 1 -2 1 0 -4 0 3; % impulse responsex = 1 -2 3 -4 3 2 1; % input sequencey = conv(h,x);n = 0:14;subplot(2,1,1);stem(n,y);xlabel(Time index n); ylabel(Amplitude);title(Output Obtained by Convolution); grid;x1 = x zeros(1,8);y1 = filter(h,1,x1);subplot(2,1,2);stem(n,y1);xlabel(Time index n); ylabel(Amplitude);title(Output Generated by Filtering); grid;8、 编写输入信号经滤波形成的系统输出信号。% Program P2_9% Generate the input sequenceclf;n = 0:299;x1 = cos(2*pi*10*n/256);x2 = cos(2*pi*100*n/256);x = x1+x2;% Compute the output sequencesnum1 = 0.5 0.27 0.77;y1 = filter(num1,1,x); % Output of System #1den2 = 1 -0.53 0.46;num2 = 0.45 0.5 0.45;y2 = filter(num2,den2,x); % Output of System #2% Plot the output sequencessubplot(2,1,1);plot(n,y1);axis(0 300 -2 2);ylabel(Amplitude);title(Output of System #1); grid;subplot(2,1,2);plot(n,y2);axis(0 300 -2 2);xlabel(Time index n); ylabel(Amplitude);title(Output of System #2); grid;9、四、本实验用到的matlab命令Stem plot sin abs cos conv clf subplot filter impz 实验二：时域连续时间信号和频域抽样理论的验证与观察 一、实验目的：1、理解并掌握信号时域采样和频率抽样理论涉及的过程和效果；2、通过编程加深理解奈奎斯特采样定理，理解不满足采样条件的对时域和频域采样造成的混叠现象。二、实验内容：1、时域的采样过程、采样定理和混叠效果；2、频域中的采样效果；3、学习buttworth模拟低通滤波器的设计命令；三、实验例程1、 连续时间信号的理想抽样及其混叠效果clf;T = 0.4;f = 25;n = (0:T:1);xs = cos(2*pi*f*n);t = linspace(-0.5,1.5,500);ya = sinc(1/T)*t(:,ones(size(n) - (1/T)*n(:,ones(size(t)*xs;plot(n,xs,o,t,ya);grid;xlabel(Time, msec);ylabel(Amplitude);title(Reconstructed continuous-time signal y_a(t);axis(0 1 -1.2 1.2)2、.频率抽样及其混叠效果clf;t = 0:0.002:50;xa = 2*t.*exp(-t);subplot(4,2,1)plot(t,xa);gridxlabel(Time, msec);ylabel(Amplitude);title(Continuous-time signal x_a(t);subplot(4,2,2)wa = 0:10/511:10;ha = freqs(2,1 2 1,wa);plot(wa/(2*pi),abs(ha);grid;xlabel(Frequency, kHz);ylabel(Amplitude);title(|X_a(jOmega)|);axis(0 5/pi 0 2);subplot(4,2,3)T = 1;n = 0:T:10;xs = 2*n.*exp(-n);k = 0:length(n)-1;stem(k,xs);grid;xlabel(Time index n);ylabel(Amplitude);title(Discrete-time signal xn);subplot(4,2,4)wd = 0:pi/255:pi;hd = freqz(xs,1,wd);plot(wd/(T*pi), T*abs(hd);grid;xlabel(Frequency, kHz);ylabel(Amplitude);title(|X(ejomega)|);axis(0 1/T 0 2)3、buttworth模拟低通滤波器的设计命令并画出该滤波器图形。clf;Fp = 3500;Fs = 4500;Wp = 2*pi*Fp; Ws = 2*pi*Fs;N, Wn = buttord(Wp, Ws, 0.5, 30,s);b,a = butter(N, Wn, s);wa = 0:(3*Ws)/511:3*Ws;h = freqs(b,a,wa);plot(wa/(2*pi), 20*log10(abs(h);gridxlabel(Frequency, Hz);ylabel(Gain, dB);title(Gain response);axis(0 3*Fs -60 5);四、本实验用到的matlab命令Butter buttord freqs freqz sinc hist实验三：离散时间信号与离散时间系统系统的频域分析 一、 实验目的：1、 掌握matlab 编写基于离散时间傅立叶变换（DTFT）、z变换和离散傅立叶变换（DFT）的程序，并通过本训练加深对这些概念和算法的理解；2、 理解DTFT、ZT和DFT的相互关系。3、 掌握离散系统的频率响应分析；4、 理解零、极点分布的概念。 5、理解帕斯瓦尔定理。二、实验内容：1、 离散时间傅立叶变换（DTFT）的定义、计算和基本性质；2、 Z变换分析；3、 系统传输函数和频率响应；4、 系统传输函数的类型；5、 系统稳定性测试；6、 离散傅立叶变换（DFT）的计算和基本性质；7、 利用FFT实现线性卷积； 8、 利用FFT显示理解帕斯瓦尔定理。三、实验例程1、用MATLAB编写离散时间信号的傅里叶变换（DTFT）并绘出8点幅频和相频曲线。% Program P3_1% Evaluation of the DTFT clf;% Compute the frequency samples of the DTFTw = -4*pi:8*pi/511:4*pi;num = 2 1;den = 1 -0.6;h = freqz(num, den, w);% Plot the DTFTsubplot(2,1,1)plot(w/pi,real(h);gridtitle(Real part of H(ejomega)xlabel(omega /pi);ylabel(Amplitude);subplot(2,1,2)plot(w/pi,imag(h);gridtitle(Imaginary part of H(ejomega)xlabel(omega /pi);ylabel(Amplitude);pausesubplot(2,1,1)plot(w/pi,abs(h);gridtitle(Magnitude Spectrum |H(ejomega)|)xlabel(omega /pi);ylabel(Amplitude);subplot(2,1,2)plot(w/pi,angle(h);gridtitle(Phase Spectrum argH(ejomega)xlabel(omega /pi);ylabel(Phase in radians);2、MATLAB验证实序列的离散时间傅立叶变换的对称关系。N = 8;k = 0:N-1;gamma = 0.5;x = exp(gamma*k);y = exp(gamma*fliplr(k);xev =0.5*(zeros(1,N-1) x+y zeros(1,N-1);xod = 0.5*(zeros(1,N-1) x-y zeros(1,N-1);X,w = freqz(x,1,512);Xev,w = freqz(xev,1,512);Xod,w = freqz(xod,1,512);Xev = exp(j*w*(N-1).*Xev;Xod = exp(j*w*(N-1).*Xod;% Verify real(X)= Xev, and imag(X)= Xodr = 0:511;w0 = -pi*r/512;X1 = freqz(x,1,w0);% Verify X = conj(X1)% real(X)= real(X1) and imag(X)= -imag(X1)% abs(X)= abs(X1) and angle(X)= -angle(X1)3、 用MATLAB计算序列Modulation Property of DTFT并给出其调制图形。% Program P3_5% Modulation Property of DTFTclf;w = -pi:2*pi/255:pi;x1 = 1 3 5 7 9 11 13 15 17;x2 = 1 -1 1 -1 1 -1 1 -1 1;y = x1.*x2;h1 = freqz(x1, 1, w);h2 = freqz(x2, 1, w);h3 = freqz(y,1,w);subplot(3,1,1)plot(w/pi,abs(h1);gridtitle(Magnitude Spectrum of First Sequence)subplot(3,1,2)plot(w/pi,abs(h2);gridtitle(Magnitude Spectrum of Second Sequence)subplot(3,1,3)plot(w/pi,abs(h3);gridtitle(Magnitude Spectrum of Product Sequence)4、用MATLAB计算序列的N点离散傅立叶变换，N取3，5，7,9。比较在时计算得到的离散时间傅立叶变换的结果。N = 16;clf;N = input(The value of N = );k = -N:N;y1=ones(1,2*N+1);y2=y1-abs(k)/N;w = 0:2*pi/255:2*pi;Y2 = freqz(y2, 1, w);Y2dft = fft(y2);k = 0:1:2*N;plot(w/pi,abs(Y2),k*2/(2*N+1),abs(Y2dft),o);xlabel(Normalized frequency);ylabel(Amplitude);4、用MATLAB研究滤波器的系统函数及其冲激响应，产生零极点图，最终显示出零点的位置。clf;b = 1 -8.5 30.5 -63;num1 = b 81 fliplr(b);num2 = b 81 81 fliplr(b);num3 = b 0 -fliplr(b);num4 = b 81 -81 -fliplr(b);n1 = 0:length(num1)-1;n2 = 0:length(num2)-1;subplot(2,2,1); stem(n1,num1);xlabel(Time index n);ylabel(Amplitude); grid;title(Type 1 FIR Filter);subplot(2,2,2); stem(n2,num2);xlabel(Time index n);ylabel(Amplitude); grid;title(Type 2 FIR Filter);subplot(2,2,3); stem(n1,num3);xlabel(Time index n);ylabel(Amplitude); grid;title(Type 3 FIR Filter);subplot(2,2,4); stem(n2,num4);xlabel(Time index n);ylabel(Amplitude); grid;title(Type 4 FIR Filter);pausesubplot(2,2,1); zplane(num1,1);title(Type 1 FIR Filter);subplot(2,2,2); zplane(num2,1);title(Type 2 FIR Filter);subplot(2,2,3); zplane(num3,1);title(Type 3 FIR Filter);subplot(2,2,4); zplane(num4,1);title(Type 4 FIR Filter);disp(Zeros of Type 1 FIR Filter are);disp(roots(num1);disp(Zeros of Type 2 FIR Filter are);disp(roots(num2);disp(Zeros of Type 3 FIR Filter are);disp(roots(num3);disp(Zeros of Type 4 FIR Filter are);disp(roots(num4);4. 用MATLAB产生序列的圆周移位,并画图观察该现象，clf;x = 0 2 4 6 8 10 12 14 16;N = length(x)-1; n = 0:N;y = circshift(x,5);XF = fft(x);YF = fft(y);subplot(2,2,1)stem(n,abs(XF);gridtitle(Magnitude of DFT of Original Sequence);subplot(2,2,2)stem(n,abs(YF);gridtitle(Magnitude of DFT of Circularly Shifted Sequence);subplot(2,2,3)stem(n,angle(XF);gridtitle(Phase of DFT of Original Sequence);subplot(2,2,4)stem(n,angle(YF);gridtitle(Phase of DFT of Circularly Shifted Sequence);5、 用MATLAB编写快速傅里叶变换，并实现偶序列的DFT的虚实部，采用图形表示。% Program P3_11% Relations between the DFTs of the Periodic Even% and Odd Parts of a Real Sequencex = 1 2 4 2 6 32 6 4 2 zeros(1,247);x1 = x(1) x(256:-1:2);xe = 0.5 *(x + x1);XF = fft(x);XEF = fft(xe);clf;k = 0:255;subplot(2,2,1);plot(k/128,real(XF); grid;ylabel(Amplitude);title(Re(DFTxn);subplot(2,2,2);plot(k/128,imag(XF); grid;ylabel(Amplitude);title(Im(DFTxn);subplot(2,2,3);plot(k/128,real(XEF); grid;xlabel(Time index n);ylabel(Amplitude);title(Re(DFTx_en);subplot(2,2,4);plot(k/128,imag(XEF); grid;xlabel(Time index n);ylabel(Amplitude);title(Im(DFTx_en);6、 编写稳定性测试程序% Program P4_4% Stability Testclf;den = input(Denominator coefficients = );ki = poly2rc(den);disp(Stability test parameters are);disp(ki);7、用MATLABD 编写基于FFT的帕斯瓦尔定理。% Program P3_12% Parsevals Relationx = (1:128) (128:-1:1);XF = fft(x);a = sum(x.*x)b = round(sum(abs(XF).2)/256)四、本实验用到的matlab命令filter Freqz impz filrfilt residuez tf2zp zp2sos zplane fft sinc zplane 实验四：IIR&FIR数字滤波器的设计一、 实验目的：1、熟悉无限冲激响应（IIR）和有限冲激响应（FIR）网络结构，对比学习模拟滤波器和数字滤波器的常用指标；2、熟悉冲激响应不变法和双线性变换法设计低通滤波器的程序编写，并深化理解；3、熟悉FIR线性相位滤波器的概念及其表示；4、熟悉FIR滤波器窗函数设计法；5、熟悉两种滤波器设计过程，并可修改其设计指标，对比输出效果。二、实验内容：1、IIR(无限冲激响应滤波器)阶数估计及其buttworth和chybyshev滤波器设计；2、冲激响应不变法和双线性变换法设计；3、FIR滤波器设计中出现频域出现吉布斯现象的来由；3、有限冲激响应滤波器窗函数设计法；三、实验计算及例程1、设低通DF的3dB带宽频率wc=0.2,止带频率ws=0.4,在 w=ws处的止带衰减20lg|H(ejws)|=-15dB,试用脉冲响应不变法（冲激不变法）设计一个Butterworth低通DF。(设采样频率fs=20kHz） Wp=input(Normalized passband edge =);Ws=input(Normalized stopband edge =);Rp=input(Passband ripple in dB =);Rs=input(Minimum stopband attenuation in dB =);N,Wn=buttord (Wp,Ws,Rp,Rs);b,a=butter(N,W n);h,omega=freqz(b,a,512);plot(omega/pi,20*log10(abs(h);grid;xlabel(omega/pi);ylabel(Gain dB);title(IIR Butterworth Lowpass Filter);2无限冲激响应滤波器阶数估计和滤波器设计 巴特沃兹带阻滤波器Ws = 0.4 0.6; Wp = 0.2 0.8; Rp = 0.4; Rs = 50;% Estimate the Filter OrderN1, Wn1 = buttord(Wp, Ws, Rp, Rs); % Design the Filternum,den = butter(N1,Wn1,stop); % Display the transfer functiondisp(Numerator Coefficients are );disp(num);disp(Denominator Coefficients are );disp(den);% Compute the gain responseg, w = gain(num,den); % Plot the gain responseplot(w/pi,g);grid axis(0 1 -60 5);xlabel(omega /pi); ylabel(Gain in dB);title(Gain Response of a Butterworth Bandstop Filter)3、吉布斯环的生成； 在FIR窗函数设计法中，频域傅立叶变换的无限振荡谐波数列叠加合成时，会产生尖锐不连续的区间，好比方波的下降及上升，这就是吉布斯现象，这种现象产生的不仅与谐波的叠加合成的数量有关，而且其变化的宽度随合成的个数增加而变窄。单位方波： 傅立叶级数表示： function y = gibbs(x,M,duty)% GIBBS% x： 待估价的值 (0x0.5)% M： 区间% duty：占空比(0-1)%if (nargin3) duty = 0.5;endp = (-M:M);for (j=1:length(x) y(j) = sum(duty*sinc(duty*p).*exp(-i*2*pi*p*abs(x(j);end正规重构：覆盖整个单位方波范围s = linspace(-0.5,0.5,200);plot(s,gibbs(s,M);4、有限冲激响应滤波器阶数估计和滤波器设计； 用凯塞窗设计一FIR低通滤波器，低通边界频率 ，阻带边界频率 ，阻带衰减 不小于50dB。 计算与例程由过渡带宽和阻带衰减 来决定凯塞窗的N和 凯塞窗设计对应的MATLAB程序为： wn=kaiser(30,4.55); nn=0:1:29; alfa=(30-1)/2; hd=sin(0.4*pi*(nn-alfa)./(pi*(nn-alfa); h=hd.*wn; h1,w1=freqz(h,1); plot(w1/pi,20*log10(abs(h1); axis(0,1,-80,10); grid;xlabel(归一化频率/p) ylabel(幅度/dB) 5、试用双线性变换法设计Butterworth低通DF。已知低通DF的3dB带宽频率，止带起始频率，在 处的止带衰减要求进行计算机MATLAB仿真并得出结果。解：使用Matlab软件仿真仿真的程序如下： clc;clear all;OmegaP=2*pi*2