江大MATLAB语音信号的采集与处理课设(4)

双线性法低通滤波器滤波后的程序：

clear all; Ft=22050; Fp=1000; Fs=1200;

[y,fs,nbits]=wavread('h15.wav'); %sound(y,fs,nbits); y=y-mean(y); as=100; ap=1;

wp=2*pi*Fp/Ft; ws=2*pi*Fs/Ft; fs1=1;

fp=2*tan(wp/2); fs=2*tan(ws/2);

[n11,wn11]=ellipord(wp,ws,ap,as,'s'); [b11,a11]=ellip(n11,ap,as,wn11,'low','s'); [num11,den11]=bilinear(b11,a11,fs1); [h,w]=freqz(num11,den11,512,Ft); subplot(2,1,1);

plot( w,20*log10(abs(h)) ); %画对数幅度谱 hold on;

xlabel( '归一化频率w/pi' ); ylabel( '幅度(dB)' ); title( 'IIR-幅度响应');

subplot(2,1,2);

plot( w,angle(h) ); hold on;

xlabel( '归一化频率w/pi' ); ylabel( '相位' ); title( 'IIR-相位响应');

f1=filter(num11,den11,y); figure(2)

subplot(2,1,1)

plot(y) %画出滤波前的时域图 title('IIR低通滤波器滤波前的时域波形'); subplot(2,1,2)

plot(f1); %画出滤波后的时域图 title('IIR低通滤波器滤波后的时域波形');

sound(f1); %播放滤波后的信号 F0=fftshift(abs(fft(f1)));

w = linspace(-Ft/2,Ft/2,length(F0)); figure(3)

y2=fftshift(abs(fft(y)));

subplot(2,1,1);

plot(w,y2); %画出滤波前的频谱图 title('IIR低通滤波器滤波前的频谱') xlabel('频率/Hz'); ylabel('幅值');

w = linspace(-Ft/2,Ft/2,length(y2)); subplot(2,1,2)

plot(w,abs(F0)); %画出滤波后的频谱图 title('IIR低通滤波器滤波后的频谱') xlabel('频率/Hz'); ylabel('幅值');

双线性法高通滤波器滤波的程序：

clear all; Ft=22050; Fp=5000; Fs=4800;

[y,fs,nbits]=wavread('h15.wav'); %sound(y,fs,nbits); y=y-mean(y); as=100; ap=1;

wp=2*pi*Fp/Ft; ws=2*pi*Fs/Ft; fs1=1;

fp=2*tan(wp/2); fs=2*tan(ws/2);

[n11,wn11]=ellipord(wp,ws,ap,as,'s'); [b11,a11]=ellip(n11,ap,as,wn11,'high','s'); [num11,den11]=bilinear(b11,a11,fs1); [h,w]=freqz(num11,den11,512,Ft); subplot(2,1,1);

plot( w,20*log10(abs(h)) ); %画对数幅度谱 hold on;

xlabel( '归一化频率w/pi' ); ylabel( '幅度(dB)' ); title( 'IIR-幅度响应');

subplot(2,1,2);

plot( w,angle(h) ); hold on;

xlabel( '归一化频率w/pi' ); ylabel( '相位' ); title( 'IIR-相位响应');

f1=filter(num11,den11,y); figure(2)

subplot(2,1,1)

plot(y) %画出滤波前的时域图 title('IIR高通滤波器滤波前的时域波形'); subplot(2,1,2)

plot(f1); %画出滤波后的时域图 title('IIR高通滤波器滤波后的时域波形');

sound(f1); %播放滤波后的信号 F0=fftshift(abs(fft(f1)));

w = linspace(-Ft/2,Ft/2,length(F0)); figure(3)

y2=fftshift(abs(fft(y))); subplot(2,1,1);

plot(w,y2); %画出滤波前的频谱图 title('IIR高通滤波器滤波前的频谱') xlabel('频率/Hz'); ylabel('幅值');

w = linspace(-Ft/2,Ft/2,length(y2)); subplot(2,1,2)

plot(w,abs(F0)); %画出滤波后的频谱图 title('IIR高通滤波器滤波后的频谱') xlabel('频率/Hz'); ylabel('幅值');

双线性法设计带通滤波器滤波的程序：

clear all; Ft=22050; Fp1=1200; Fp2=3000; Fs1=1000; Fs2=3200;

[y,fs,nbits]=wavread('h15.wav'); %sound(y,fs,nbits); y=y-mean(y); Fp=[Fp1,Fp2]; Fs=[Fs1,Fs2]; as=100; ap=1; T=2;

wp1=2*pi*Fp1/Ft; wp2=2*pi*Fp2/Ft; ws1=2*pi*Fs1/Ft; ws2=2*pi*Fs2/Ft;

Wp1=(2/T)*tan(wp1/2); Wp2=(2/T)*tan(wp2/2);

Ws1=(2/T)*tan(ws1/2); Ws2=(2/T)*tan(ws1/2)

W0=Wp1*Wp2; w0=sqrt(W0);

BW=Wp2-Wp1; %带通滤波器的通带宽度 lp=1; %归一化处理

ls=Ws1*BW/(W0-Ws1^2); [N,Wn]=ellipord(lp,ls,ap,as,'s'); [B,A]=ellip(N,1,40,Wn,'s'); [BT,AT]=lp2bp(B,A,w0,BW); [num,den]=bilinear(BT,AT,0.5); [z,p,k]=tf2zp(num,den);

[h,w]=freqz(num,den,512,Ft); figure(1)

subplot(2,1,1);

plot( w,20*log10(abs(h)) ); %画对数幅度谱 hold on;

xlabel( '频率w/pi' ); ylabel( '幅度(dB)' ); title( 'IIR-幅度响应'); subplot(2,1,2);

plot( w,angle(h) ); hold on;

xlabel( '频率w/pi' ); ylabel( '相位' ); title( 'IIR-相位响应');

f1=filter(num,den,y); figure(2)

subplot(2,1,1)

plot(y) %画出滤波前的时域图 title('IIR带通滤波器滤波前的时域波形'); subplot(2,1,2)

plot(f1); %画出滤波后的时域图 title('IIR带通滤波器滤波后的时域波形');

sound(f1); %播放滤波后的信号 F0=fftshift(abs(fft(f1)));

w = linspace(-Ft/2,Ft/2,length(F0)); figure(3)

y2=fftshift(abs(fft(y))); subplot(2,1,1);

plot(w,y2); %画出滤波前的频谱图 title('IIR带通滤波器滤波前的频谱') xlabel('频率/Hz'); ylabel('幅值');