数字信号matalab实验.doc

实验二 IIR滤波器设计一、 实验目的1、 掌握冲激响应法和双线性变换法设计IIR滤波器的原理及具体设计方法，熟悉用双线性设计法设计低通、带通和高通IIR数字滤波器的计算机程序。2、 熟悉模拟Butterworth滤波器的设计，掌握冲激响应法和双线性变换法设计数字IIR滤波器的方法。二、 实验内容1、 不同阶次模拟巴特沃兹滤波器的频率响应如下： 由图可知，不同阶次的所对应的滤波器的幅度走向基本一样，但是过渡带却有着明显的不同，阶次越高，滤波器的过渡带越窄。2、根据模拟滤波器指标，设计低通巴特沃兹滤波器设计指标为：通带截止频率fp = 5kHz, 通带最大衰减 p =2dB，阻带截止频率 fs= 12kHz, 阻带最小衰减s =30dB。设计程序如下： clear all; fp=5e3; ap=2; fs=12e3; as=30; %*求阶次N*% k_sp=sqrt(10(as/10) - 1)/sqrt(10(ap/10) - 1); a_sp = 2*pi*fs/(2*pi*fp); N=ceil( log(k_sp)/log(a_sp) ); p=zeros(1,N); for i=1:N;k = i-1;p(i) = exp( j*pi/2 + j*pi*(2*k+1)/2/N ) ;end z=; b1 a1 = zp2tf( z, p, 1 ); fp1=fp* ( 10(ap/10)-1 )(-1/2/N) ; fs1=fs* ( 10(as/10)-1 )(-1/2/N); figure(1); x= 0:50:1.5*fs ; H1 =freqs(b1, a1, x/fp1); H2 =freqs(b1, a1,x/fs1); plot( x, 20*log10(abs(H1); xlabel( 频率/Hz ); ylabel( 振幅/db ); title( 巴特沃兹滤波器对数幅度谱 N=,num2str(N) ); figure(2); x1= 0: 0.01 :2 ; H3 =freqs(b1, a1, x1); plot( x1 , abs(H3) ); xlabel(归一化频率); ylabel(幅度); title( 归一化之后的低通巴特沃兹滤波器的频率响应: );归一化之前低通巴特沃兹滤波器的频率响应:归一化之后的低通巴特沃兹滤波器的频率响应为: 3、已知模拟滤波器传输函数H(s)为按照不同采样频率Fs1=1Hz, Fs2=10Hz 使用脉冲响应不变法将H(s)转换为H(z)数字IIR滤波器。不同采样采样频率对应的IIR滤波器的频率响应见下：程序如下： clear all; b = 0 0 0.5012; a = 1 0.6449 0.7079; f= 0:0.5:30; H = freqs(b, a, f); figure(1); plot( f,20*log10(abs(H) ); xlabel( 角频率( rad/s ) ); ylabel( 振幅/dB ); title(模拟滤波器频率响应); Fs1 = 1; Fs2 = 10; b1, a1 = impinvar( b, a, Fs1 ); b2, a2 = impinvar( b, a, Fs2 ); f1= 0:0.1:2*pi; H1 = freqz( b1, a1, f1 ); H2 = freqz( b2, a2, f1 ); figure(2); plot( f1/pi, 20*log10(abs(H1) ); xlabel( 角频率/Hz ); ylabel( 幅频响应/dB) ); title( 不同采样频率对于数字IIR频响的影响 ); hold on; plot( f1/pi, 20*log10(abs(H2),r ); legend( Fs1= num2str(Fs1) Hz, Fs2= num2str(Fs2) Hz );从图可知由不同的采样频率设计出的IIR滤波器的频率响应明显是不一样的，采样频率越高，所设计出的频率计在pi处衰减的越大，滤波器越接近模拟滤波器。4 给定数字IIR指标，用脉冲响应不变法设计IIR。 在不同的Fs下，所得到的模拟滤波器和数字滤波器的频率响应见下：程序如下：clear all;freq_p = 0.2*pi; freq_s = 0.3*pi;alpha_p = 1;alpha_s = 15; Fs_array = 1 10 ;Color_array = b,k ;for i=1:2 color_used = Color_array(i); Fs = Fs_array(i); T = 1/Fs;analog_freq_p = freq_p/T %*数字频率w到模拟角频率W的转换 W = w/T% analog_freq_s = freq_s/T %*阶次N*% k_sp = sqrt(10(alpha_p/10) - 1)/sqrt(10(alpha_s/10) - 1); lamda_sp = 2*pi*analog_freq_s/(2*pi*analog_freq_p); N = ceil( -log(k_sp)/log(lamda_sp) ) Wc = analog_freq_p* ( 10(analog_freq_p/10)-1 )(-1/2/N) %3dB截至频率 %*求极点*% p = zeros(1,N); for i=1:N k = i-1; p(i) = exp( j*pi/2 + j*pi*(2*k+1)/2/N ) ; endz = ; b a = zp2tf( z, p, 1 ); % 根据p=s/Wc 把归一化传输函数H(p)的系数转换为实际的模拟传输函数H(s) bs = zeros(1,length(b); as = zeros(1,length(a); temp =N:-1:0; temp = Wc . temp; bs = b ./ temp; as = a ./ temp; temp = as(1); as = as/temp ;%验证与书上结果是否一致 bs = bs/temp; figure(1); freq_analog = linspace( 0, 2*analog_freq_s ); H_analog = freqs(bs,as,freq_analog); plot( freq_analog, 20*log10(abs(H_analog) ,color_used); hold on; xlabel( 模拟角频率(rad/s) ); ylabel(振幅(dB); title( 模拟滤波器的频率响应 ); % 根据H(s)获得H（z），使用脉冲响应不变法 bz, az = impinvar( bs, as, Fs ); %根据模拟滤波器系数b,a，按照采样率fs转换为系数为bz,az的数字IIR滤波器 freq_digital = 0:0.1:2*pi-0.1; H = freqz( bz, az, freq_digital ); figure(2); plot( freq_digital/pi, 20*log10(abs(H) ,color_used ); grid on; hold on; xlabel( 数字频率/pi); ylabel( 振幅(dB); title( 数字IIR滤波器频率响应 ); a = 1;end% 当采样频率变为2Hz时，求得出的数字IIR频率响应figure(1),legend(Fs=num2str(Fs_array(1)Hz,Fs= num2str(Fs_array(2) Hz );figure(2),legend(Fs=num2str(Fs_array(1)Hz,Fs= num2str(Fs_array(2) Hz ); 由上可知在不同的采样频率下，由数字指标所得到模拟指标不一样，因此导致所设计的模拟滤波器的频率响应也不一样。但是，从结果看由模拟滤波器转化为的数字滤波器差别不大，符合所给的数字指标。5、 用双线性变换法设计IIR，并同脉冲响应不变法比较设计的结果。数字滤波器指标要求如下：设计程序如下： clear all; fp = 0.2*pi; ap = 1; fs = 0.3*pi; as = 15; % *采样频率*% Fs = 1; T = 1/Fs; analog_fp = 2/T * tan(fp/2); analog_fs = 2/T * tan(fs/2); %* 求节数N*% a_sp = 2*pi*analog_fs/(2*pi*analog_fp); k_sp = sqrt(10(as/10) - 1)/sqrt(10(ap/10) - 1); N = ceil( log(k_sp)/log(a_sp) ); Wc = analog_fp* ( 10(analog_fp/10)-1 )(-1/2/N) ; %3dB截至频率 p = zeros(1,N); for i=1:N k = i-1; p(i) = exp( j*pi/2 + j*pi*(2*k+1)/2/N ) ; end z = ; b a = zp2tf( z, p, 1 ); %得到传递函数 bs = zeros(1,length(b); as = zeros(1,length(a); temp =N:-1:0; temp = Wc . temp; bs = b ./ temp; as = a ./ temp; temp = as(1); as = as/temp; bs = bs/temp; figure(1); freq_analog = linspace( 0, 2*analog_fs ); H_analog = freqs(bs,as,freq_analog); plot( freq_analog, 20*log10(abs(H_analog) ,k); grid on; xlabel( 模拟角频率(rad/s) ); ylabel(振幅(dB); title( 模拟滤波器的频率响应 ); % 根据模拟滤波器的传输函数得到数字滤波器的传输函数 bz az = bilinear( bs, as, Fs ); freq_digital = 0:0.1:2*pi-0.1; H = freqz( bz, az, freq_digital ); figure(2); plot( freq_digital/pi, 20*log10(abs(H) ); grid on; xlabel( 数字频率/pi); ylabel( 振幅(dB); title( 数字IIR滤波器频率响应 ); 与题目给的要求对比，可以得之所设计的滤波器符合题目的要求。三、 实验结论由冲激响应不变法所设计的滤波器在直流附近的数字域与模拟域的幅频特性相近，但双线性变换法是以严重的失真换取了频率的不混叠，这一般用于具有片段常数特性的滤波器。实验三 FIR滤波器设计一、 实验目的1、 熟悉FIR滤波器设计的基本方法；2、 掌握窗函数法设计FIR滤波器的原理和方法；3、 熟悉线性相位FIR滤波器的幅频特性和相位特性；4、 了解不同窗函数对滤波器性能的响应。二、 实验内容1、 验证窗函数N变化时，验证其频谱主瓣副瓣比、主瓣宽度的变化。以矩形窗为例来说明N为20 时：程序如下： clear all N=20; n=ones(1,N); f1=fft(n,128); f1=f1/max(f1); x=(0:1:127)/128*2*pi; figure(1); plot(x,20*log10(abs(f1); xlabel(频率); ylabel(振幅/db); title(N=20);图形如下： N为30 时：程序如下： clear all N=30; n=ones(1,N); f1=fft(n,128); f1=f1/max(f1); x=(0:1:127)/128*2*pi; figure(1); plot(x,20*log10(abs(f1); xlabel(频率); ylabel(振幅/db); title(N=30); 从图中可以算出，当N=20时主瓣副瓣幅度比为4.5704，主瓣宽为0.63pi，约为4*pi/N; 当N=30时主瓣副瓣幅度比为4.59，主瓣宽为0.43pi，约为4*pi/N;故可得：矩形窗函数的N增大时，验证其其频谱主瓣副瓣幅度比基本不会发生变化，而主瓣宽度将会变窄。这说明，当用矩形窗函数设计滤波器时，增大N不能使得阻带衰减增大，但能够减小过渡带。2、用窗函数法设计线性相位FIR低通通带截止频率wp=0.5PI, 阻带截止频率ws=0.6PI, 阻带衰减不小于40dB,通带衰减不大于3dB.a) 选取Hanning和Hamming就可以满足设计的要求。b) 根据过渡带的要求选择窗的长度。过渡带是0.1 PI,海明和汉宁窗的主瓣宽度均是8 PI/N的，从而选择的窗的长度都是大于等于80的，选择N大于等于80就可以了，在这里选取N=81.用Hanning和Hamming窗设计的滤波器见下面：程序如下： clear all wp = 0.5*pi; ws = 0.6*pi; width = ws - wp; % *Hanning窗*% N = ceil( 8*pi/width );%* 截止频率*% Wn = (wp + ws)/2; if mod(N,2)=0N = N + 1;end window_hanning = hanning(N); b_hanning = fir1( N-1, Wn/pi, window_hanning ); norm = 0:511/512; H_hanning = freqz( b_hanning, 1, 512); subplot(2,1,1); plot(norm,20*log10(abs(H_hanning) ); hold on; xlabel( 归一化频率w/pi ); ylabel( 振幅/dB ); title( 幅频响应); subplot(2,1,2); plot( norm,angle(H_hanning) ); hold on; xlabel( 归一化频率w/pi ); ylabel( 相位 ); title( 相频响应); %* Hamming窗*% N1= ceil( 8*pi/width); Wn = (wp + ws)/2; %截止频率 if mod(N1,2)=0 N1= N1 + 1;end window_hamming = hamming(N1); b_hamming = fir1( N1-1, Wn/pi, window_hamming ); H_hamming = freqz( b_hamming, 1, 512); subplot(2,1,1); plot(norm, 20*log10(abs(H_hamming),r ); subplot(2,1,2); plot( norm,angle(H_hamming),r ); legend( Hanning num2str(N) 阶 ,Hamming num2str(N1) 阶 );图如下： 很明显从图中可以看出汉宁窗和海明窗都符合要求，此时窗长必须为奇数。3、 用窗函数法设计线性相位FIR高通，通带截止频率wp=0.8PI, 阻带截止频率ws=0.7PI, 阻带衰减不小于30dB,通带衰减不大于3dB. 1）根据阻带衰减的要求选择窗函数，可以选择汉宁窗（最大衰减是44dB）2)根据过渡带的要求选择窗的长度。过渡带是0.1 PI,海明和汉宁窗的主瓣宽度均是8 PI/N的，从而选择的窗的长度都是大于等于80的，选择N=81（必须是奇数）;程序如下： clear all; wp = 0.8*pi; ws = 0.7*pi; width = wp - ws;%* Hanning窗*% N = ceil( 8*pi/width ); Wn = (wp + ws)/2; if mod(N,2)=0 N = N + 1;end window_hanning = hanning(N); b_hanning = fir1( N-1, Wn/pi, high,window_hanning ); norm = 0:511/512; H_hanning = freqz( b_hanning, 1, 512); subplot(2,1,1); plot( norm,20*log10(abs(H_hanning) ); hold on; xlabel( 归一化频率w/pi ); ylabel( 幅度(