ch3_9_42.m

用MATLAB作的一个滤波程序

clear all;
figure(2)
[y1,Fs,bits]=wavread('Windows');%采点106150
%抽样频率Fs=22050
a1=y1(:,[1])
b1=a1';%通道1采点的行向量
c1=b1(1:3:end);%3倍抽取 35384点

a2=y1(:,[2])
b2=a2';%通道1采点的行向量
c2=b2(1:3:end);%3倍抽取  35384点
sound([c1',c2'])
pause;
%设计双带阻滤波器滤除500HZ和1000HZ的频谱
k=0:35383;
%c1=c1(1:3000);
%c2=c2(1:3000);
d1=c1+0.5*sin(2*pi*500*k./22050)+0.5*sin(2*pi*1000*k./22050);
d2=c2+0.5*sin(2*pi*500*k./22050)+0.5*sin(2*pi*1000*k./22050);
sound([d1',d2']);%加噪声后试听

%窗口法设计双带通滤波器
figure(2);


f=[490,495,505,510,980,990,1010,1020];
f1=f./Fs.*2*pi
a=[1,0,1,0,1];
Ap=1;As=80;
devp = (10^(Ap/20)-1)/(10^(Ap/20)+1); 
devr=10^(-As/20);
devs=[devp,devr,devp,devr,devp];
[M,Wc,beta,ftype]=kaiserord(f,a,devs,Fs);
M=mod(M,2)+M;
h=fir1(M,Wc,ftype,kaiser(M+1,beta))
;
w=linspace(0,pi,512);
mag=freqz(h,1,w);
subplot(221)
plot(w/pi,20*log10(abs(mag)));
xlabel('Frequency ');
ylabel('Magnitude (dB)');
title('窗口法设计双带阻滤波器');
axis([0 0.2 -80 1]);
grid;

%用滤波器进行滤波
%通道1
z1=filter(h,1,d1);


%带阻滤波后信号的频谱图

fl=fft(z1,1024);
L=0:1023;
p=length(L);
X=fftshift(fl);
w1=-11025+(0:p-1)*2*11025./p;
subplot(222);
plot(w1,abs(X));
xlabel(' frequency  HZ');
ylabel('magnitude  ');
axis([0 3000 0 0.01]);
title('通道1带阻滤波后用FFT得到幅频响应')
grid;

%通道2
z2=filter(h,1,d2);


%带阻滤波后信号的频谱图

fl=fft(z2,1024);
L=0:1023;
p=length(L);
X=fftshift(fl);
w1=-11025+(0:p-1)*2*11025./p;
subplot(223);
plot(w1,abs(X));
xlabel(' frequency  HZ');
ylabel('magnitude  ');
axis([0 3000 0 0.01]);
title('通道2带阻滤波后用FFT得到幅频响应')
grid;

%试听滤波后的声音
sound([z1',z2']);