ch3_9_4.m

用MATLAB作的一个滤波程序

clear all;
figure(1);
%读入信号
[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']);%重新组合后试听
  
%分析c1,c2的频谱

f=fft(c1,1024);
L=0:1023;
p=length(L);
X=fftshift(f);
w1=-11025+(0:p-1)*2*11025./p;
subplot(221);
plot(w1,abs(X));
xlabel(' frequency  HZ');
ylabel('magnitude  ');
axis([-6000 6000 0 25]);
title('用FFT得到通道1的幅频响应')
grid;

f=fft(c2,1024);
L=0:1023;
p=length(L);
X=fftshift(f);
w1=-11025+(0:p-1)*2*11025./p;
subplot(222);
plot(w1,abs(X));
xlabel(' frequency  HZ');
ylabel('magnitude  ');
axis([-6000 6000 0 20]);
title('用FFT得到通道2的幅频响应')
grid;

pause
%加噪声
k=0:35383;
d1=c1+0.1*sin(2*pi*500*k./22050)+0.1*sin(2*pi*1000*k./22050);
d2=c2+0.1*sin(2*pi*500*k./22050)+0.1*sin(2*pi*1000*k./22050);
sound([d1',d2']);%加噪声后试听

%分析d1,d2的频谱

f=fft(d1,1024);
L=0:1023;
p=length(L);
X=fftshift(f);
w1=-11025+(0:p-1)*2*11025./p;
subplot(223);
plot(w1,abs(X));
xlabel(' frequency  HZ');
ylabel('magnitude  ');
axis([0 4000 0 25]);
title('用FFT得到通道1加噪的幅频响应')
grid;

f=fft(d2,1024);
L=0:1023;
p=length(L);
X=fftshift(f);
w1=-11025+(0:p-1)*2*11025./p;
subplot(224);
plot(w1,abs(X));
xlabel(' frequency  HZ');
ylabel('magnitude  ');
axis([0 4000 0 20]);
title('用FFT得到通道2的加噪幅频响应')
grid;