matlab 数字信号处理应用实例,数字信号处理Matlab_实现实例

数字信号处理Matlab_实现实例

数字信号处理 Matlab 实现实例第１章 离散时间信号与系统 例 1-1 用 MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的离散卷积。 解 MATLAB 程序如下： a=[-2 0 1 -1 3]; b=[1 2 0 -1]; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel( n ); ylabel( 幅度 ); 图 1.1给出了卷积结果的图形，求得的结果存放在数组 c中为：{-2 -4 1 3 1 5 1 -3}。 例 1-2 用 MATLAB计算差分方程 当输入序列为 时的输出结果 。 解 MATLAB 程序如下： N=41;a=[0.8 -0.44 0.36 0.22];b=[1 0.7 -0.45 -0.6];x=[1 zeros(1,N-1)];k=0:1:N-1;y=filter(a,b,x);stem(k,y)xlabel( n );ylabel( 幅度 )图 1.2 给出了该差分方程的前 41个样点的输出，即该系统的单位脉冲响应。 例 1-3 用 MATLAB计算例 1-2差分方程 所对应的系统函数的 DTFT。 解 例 1-2差分方程所对应的系统函数为： 1230.840.36.()75zzHz其 DTFT为230.840.36.0()175jjjj eee用 MATLAB计算的程序如下：k=256;num=[0.8 -0.44 0.36 0.02];den=[1 0.7 -0.45 -0.6];w=0:pi/k:pi;h=freqz(num,den,w);subplot(2,2,1);plot(w/pi,real(h));gridtitle( 实部 )xlabel( \omega/\pi );ylabel( 幅度 )subplot(2,2,2);plot(w/pi,imag(h));gridtitle( 虚部 )xlabel( \omega/\pi );ylabel( Amplitude )subplot(2,2,3);plot(w/pi,abs(h));gridtitle( 幅度谱 )xlabel( \omega/\pi );ylabel( 幅值 )subplot(2,2,4);plot(w/pi,angle(h));gridtitle( 相位谱 )xlabel( \omega/\pi );ylabel( 弧度 ) 第２章 离散傅里叶变换及其快速算法例 2-1 对连续的单一频率周期信号 按采样频率 采样，截取长度 N分别选 N =20和 N =16，观察其 DFT结果的幅度谱。 解 此时离散序列 ，即 k=8。用 MATLAB计算并作图，函数 fft用于计算离散傅里叶变换 DFT，程序如下： k=8; n1=[0:1:19]; xa1=sin(2*pi*n1/k); subplot(2,2,1) plot(n1,xa1) xlabel( t/T );ylabel( x(n) ); xk1=fft(xa1);xk1=abs(xk1); subplot(2,2,2) stem(n1,xk1) xlabel( k );ylabel( X(k) ); n2=[0:1:15]; xa2=sin(2*pi*n2/k); subplot(2,2,3) plot(n2,xa2) xlabel( t/T );ylabel( x(n) ); xk2=fft(xa2);xk2=abs(xk2); subplot(2,2,4) stem(n2,xk2) xlabel( k );ylabel( X(k) );计算结果示于图 2.1，(a)和(b)分别是 N=20时的截取信号和 DFT结果，由于截取了两个半周期，频谱出现泄漏；(c) 和(d) 分别是 N=16时的截取信号和 DFT结果，由于截取了两个整周期，得到单一谱线的频谱。上述频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏。 例 2-2 用 FFT计算两个序列 的互相关函数 。 解 用 MATLAB计算程序如下： x=[1 3 -1 1 2 3 3 1]; y=[2 1 -1 1 2 0 -1 3]; k=length(x); xk=fft(x,2*k);yk=fft(y,2*k); rm=real(ifft(conj(xk).*yk)); rm=[rm(k+2:2*k) rm(1:k)]; m=(-k+1):(k-1); stem(m,rm) xlabel( m ); ylabel( 幅度 );其计算结果如图 2.2所示。例 2-3计算两个序列的的互相关函数，其中 x(n)={2 3 5 2 1 –1 0 0 12 3 5 3 0 –1 –2 0 1 2} y(n)=x(n-4)+e(n), e(n)为一随机噪声，在 MATLAB中可以用随机函数 rand产生 解 用 MATLAB计算程序如下： x=[2 3 5 2 1 -1 0 0 12 3 5 3 0 -1 -2 0 1 2]; y=[0 0 0 0 2 3 5 2 1 -1 0 0 12 3 5 3 0 -1 -2 0 1 2]; k=length(y); e=rand(1,k)-0.5; y=y+e; xk=fft(x,2*k); yk=fft(y,2*k); rm=real(ifft(conj(xk).*yk)); rm=[rm(k+2:2*k) rm(1:k)]; m=(-k+1):(k-1); stem(m,rm)xlabel( m ); ylabel( 幅度 ); 计算结果如图 2.3(a)，我们看到最大值出现在 m=4处，正好是 y(n)对于 x(n)的延迟。2. 3(b)是 x(n)的自相关函数，他和 y(n)的区别除时间位置外，形状也略不同，这是由于 y(n)受到噪声的干扰。第 3章 无限长单位脉冲响应(IIR)滤波器的设计方法例 3-1 设采样周期 T=250μs(采样频率 fs =4kHz)，用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其 3dB边界频率为 fc =1kHz。 [B,A]=butter(3,2*pi*1000, s ); [num1,den1]=impinvar(B,A,4000); [h1,w]=freqz(num1,den1); [B,A]=butter(3,2/0.00025, s ); [num2,den2]=bilinear(B,A,4000); [h2,w]=freqz(num2,den2); f=w/pi*2000; plot(f,abs(h1), -. ,f,abs(h2), - ); grid; xlabel( 频率/Hz ) ylabel( 幅值/dB ) 程序中第一个 butter的边界频率 2π×1000，为脉冲响应不变法原型低通滤波器的边界频率；第二个 butter的边界频率 2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图 3.1给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。脉冲响应不变法由于混叠效应，使