MATLAB处理语音信号滤波程序.doc

MATLAB处理语音信号滤波程序

巴特沃思滤波器滤波
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 参数初始化 %%%%%%%%%%%%%%%%
clear all; % 关闭以往仿真的数据栈
A=5; %噪声的振幅
u1=0.01;%噪声振幅的衰减系数
u2=0.02;
u3=0.03;
f1=50;%多频噪声频率
f2=200;
f3=800;
[Y,fs,bits]=wavread('0601331007.wav') %读出信号1，采样率和采样位数
% sound(Y,fs);
% [X,Fs,bits]=wavread('0601331007.wav'); %读出信号2，采样率和采样位数
n=length(Y); %求采样信号1的长度
% m=length(X); %求采样信号2的长度
t=0:1/fs:(n-1)/fs; %求采样信号1时域上的采样点数
% t1=0:1/Fs:(m-1)/fs; %求采样信号2时域上的采样点数
% % % % % % t=(0:n-1)/fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 构造噪声 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y1=u1*A*sin(2*pi*f1*t)+u1*A*sin(2*pi*f2*t)%两两噪声相加
y2=u1*A*sin(2*pi*f3*t)
y=y1+y2
% % % % % % Z1=u1*A*sin(2*pi*f1*t1)+u1*A*sin(2*pi*f2*t1)%两两噪声相加
% % % % % % Z2=u1*A*sin(2*pi*f3*t1)
% % % % % % Z=Z1+Z2
% figure(1) %画出噪声的时域图，取200个点
% % % % % % % % % % plot(t ,y);grid on;
% plot(t(1:200),y(1:200));grid on;
%
% xlabel('时间 (t)');
% ylabel('幅度（y)');
% title('噪声的时域图');
%
% Fy=fft(y,n); %对噪声进行傅立叶变换
% Fy1=abs(Fy);
% n1=floor(n/2); %对zhang求采样点数的一半
% n2=floor(m/2); %对zhanglow求采样点数的一半
% f=(0:n1)*fs/n; %时域上的采样点数
% f1=(0:n2)*Fs/m; %时域上的采样点数
%
% figure(2) %画出噪声的频谱图
% plot(f,Fy1(1:n1+1));grid on;
% xlabel('频率（f）');
% ylabel(' 幅度（Fy1）');
% title('噪声的频谱图');
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 加噪声前后的时域图比较 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Y=Y(:,1);%此处假设声音是双声道的，只取单声道作分析
% X=X(:,1);%此处假设声音是双声道的，只取单声道作分析
Y1=y+Y'; %对采样信号进行转置再加上噪声
% pause(10); % 暂停10s
sound(Y1,fs); % 读出加噪声后的采样信号
% figure(3) %画出加噪声前后的时域比较图
% subplot(2,1,1);plot(t,Y);grid on;
% plot(t(1:1000),Y(1:1000));grid on;
% xlabel('时间 (t)');
% ylabel('幅度（Y)');
% title('加噪声前的时域图');
%
% subplot(2,1,2);plot(t,Y1);grid on;
% % % % % % % % % % plot(t(1:1000),Y1(1:1000));grid on;
% xlabel('时间 (t)');
% ylabel('幅度（Y1)');
% title('加噪声后的时域图');
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 加噪声前后的频域图比较 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% FY=fft(Y1,n); %对加噪声后的采样信号1进行傅立叶变换
% FY1=abs(FY);
% FY2=fft(Y,n); %对加噪声前的采样信号1进行傅立叶变换
% FY3=abs(FY2);
% FX=fft(X,m); %对加噪声前的采样信号2进行傅立叶变换
% FX4=abs(FX);
%
% figure(4) %画出加噪声前后的频域比较图
% subplot(2,1,1);
% plot(f,FY3(1:n1+1));grid on;
% xlabel('频率（f）');
% ylabel(' 幅度（FY3）');
% title('加噪声前的频谱图');
%
% subplot(2,1,2);
% plot(f,FY1(1:n1+1));grid on;
% xlabel('频率（f）');
% ylabel(' 幅度（FY1）');
% title('加噪声后的频谱图');
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 不同音调相同音量信号的频域、时域图比较 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% figure(5)
% subplot(2,2,1);plot(f,FY3(1:n1+1));grid on;
% xlabel('频率（f）');ylabel(' 幅度（FY3）');title('zhang加噪声前的频谱图');
% subplot(2,2,3);plot(f1,FX4(1:n2+1));grid on;
% xlabel('频率（f）');ylabel(' 幅度（FX4）');title('zhanglow加噪声前的频谱图');
% subplot(2,2,2);plot(t,Y);grid on;
% xlabel('时间 (t)');ylabel('幅度（Y)');title('zhang加噪声前的时域图');
% subplot(2,2,4);plot(t1,X);grid on;
% xlabel('时间 (t1)');ylabel('幅度（X)');title('zhanglow加噪声前的时域图');
% % sound(Y1,fs);%读出加噪后的信号

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 构造巴特沃思滤波器 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(6); %画出zhang加噪声后的时域图

plot(t,Y1);
xlabel('时间 (t)');ylabel('幅度（Y1)');title('lishuai加噪声后的时域图');grid on;
% axis([0.04 0.06 -3 3]);

%%%%%%%%%%%%%% 构造800Hz的带阻滤波器 %%%%%%%%%%%%%

f0_stop1=800; %
fc=20; % 设定过渡带的频率
Rp=1; % 通带区的波纹系数
Rs=30; % 阻带区的波纹系数
wp_stop1=[(f0_stop1-fc)/(fs/2) (f0_stop1+fc)/(fs/2)]; % 通带的拐角频率
ws_stop1=[(f0_stop1-5*fc)/(fs/2) (f0_stop1+5*fc)/(fs/2)]; % 阻带的拐角频率
[N_stop1,wc_stop1]=buttord(wp_stop1,ws_stop1,Rp,Rs,'s'); % 求出巴特沃思滤波器的阶数N及频率参数wc
[num_stop1,den_stop1]=butter(N_stop1,wc_stop1,'stop'); % 求出巴特沃斯带阻数字滤波器的传递函数模型系数
[h_stop1,w_stop1]=freqz(num_stop1,den_stop1,fs) % 求出离散系统频率响应的数值
figure(7);
subplot(2,3,1); % 画出带阻滤波器的幅频特性图
plot(w_stop1*fs/(2*pi),20*log10(abs(h_stop1)));
xlabel('频率（w_stop1）');ylabel('幅度（h_stop1)');title('带阻滤波器在800Hz处的幅频特性图');grid on;
B_stop1=filter(num_stop1,den_stop1,Y1); %对含噪信号Y1进行带阻滤波
% pause(10);
% sound(B_stop1,fs); % 读出滤去800Hz噪声后的采样信号
subplot(2,3,4); % 画出带阻滤波器滤去800Hz噪声后的时域图
plot(t,B_stop1);
xlabel('时间 (t)）');ylabel('幅度（B_stop1)');title('用带阻滤波器滤去800Hz噪声后的时域图');grid on;
% axis([0.04 0.06 -30 3]);
%%%%%%%%%%%%%% 构造200Hz的带阻滤波器 %%%%%%%%%%%%%%%%%%%%%%
f0_stop2=200;
wp_stop2=[(f0_stop2-fc)/(fs/2) (f0_stop2+fc)/(fs/2)]; % 通带的拐角频率
ws_stop2=[(f0_stop2-5*fc)/(fs/2) (f0_stop2+5*fc)/(fs/2)]; % 阻带的拐角频率
[N_stop2,wc_stop2]=buttord(wp_stop2,ws_stop2,Rp,Rs,'s'); % 求出巴特沃思滤波器的阶数N及频率参数wc
[num_stop2,den_stop2]=butter(N_stop2,wc_stop2,'stop'); % 求出巴特沃斯带阻数字滤波器的传递函数模型系数
[h_stop2,w_stop2]=freqz(num_stop2,den_stop2,fs) % 求出离散系统频率响应的数值
subplot(2,3,2);
plot(w_stop2*fs/(2*pi),20*log10(abs(h_stop2)));
xlabel('频率（w_stop2）');ylabel('幅度（h_stop2)');title('带阻滤波器在200Hz处的幅频特性图');grid on;
B_stop2=filter(num_stop2,den_stop2,B_stop1);%对含噪信号Y1进行带阻滤波
% pause(10);
% sound(B_stop2,fs); % 读出滤去800Hz和200Hz噪声后的采样信号
subplot(2,3,5); % 画出带阻滤波器滤去800Hz和200Hz噪声后的时域图
plot(t,B_stop2);
xlabel('时间 (t)）');ylabel('幅度（B_stop2)');title('带阻滤波器滤去800Hz和200Hz噪声后的时域图');grid on;
% axis([0.04 0.06 -3 3]);

%%%%%%%%%%%%%% 构造50Hz的带阻滤波器 %%%%%%%%%%%%%%%%%%%%%%
f0_stop3=50
ws_stop3=(f0_stop3-fc)/(fs/2) % 阻带的拐角频率
wp_stop3=(f0_stop3+5*fc)/(fs/2) % 通带的拐角频率
Wws=0.0075;
Wwp=0.0375;

[N_stop3,wc_stop3]=buttord(Wwp,Wws,Rp,Rs,'s'); % 求出巴特沃思滤波器的阶数N及频率参数wc
[num_stop3,dem_stop3]=butter(N_stop3,wc_stop3,'high'); % 求出巴特沃斯高通数字滤波器的传递函数模型系数
[h_stop3,w_stop3]=freqz(num_stop3,dem_stop3,fs); % 求出离散系统频率响应的数值
subplot(2,3,3); %画出高通滤波器的幅频特性图
plot(w_stop3*fs/(2*pi),20*log10(abs(h_stop3)));
xlabel('频率（w_stop3）');ylabel('幅度（h_stop3)');title('高通滤波器在50Hz处的幅频特性图');grid on;
B_stop3=filter(num_stop3,dem_stop3,B_stop2);
subplot(2,3,6); % 画出高通滤波器滤去800Hz和200Hz和50Hz噪声后的时域图
plot(t,B_stop3);
xlabel('时间 (t)）');ylabel('幅度（B_stop3)');title('滤去800Hz和200Hz和50Hz噪声后的时域图');grid on;
% axis([0 10 -60000 60000]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%用不同的滤波器滤去相噪声频率后的时域比较图 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(8);
subplot(2,2,1);plot(t,Y1);grid on;
xlabel('时间 (t)');ylabel('幅度（Y1)');title('加噪声后的时域图');
subplot(2,2,2);plot(t,B_stop1);
xlabel(' 时间 (t)）');ylabel('幅度（B_stop1)');
title('用带阻滤波器滤去800Hz噪声后的时域图');grid on;
subplot(2,2,3);plot(t,B_stop2);
xlabel(' 时间 (t)）');ylabel('幅度（B_stop2)');
title( ' 带阻滤波器滤去800Hz和200Hz噪声后的时域图');grid on;
subplot(2,2,4);plot(t,B_stop3);
xlabel(' 时间 (t)）');ylabel('幅度（B_stop3)');
title( '高通滤波器滤去800Hz和200Hz和50Hz噪声后的时域图');grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%% 用不同的滤波器滤去相应噪声频率后的频域比较 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(9);

FY=fft(Y1,n); %对加噪声后的采样信号1进行傅立叶变换
FY1=abs(FY);
n1=floor(n/2);
f=(0:n1)*fs/n;
FY_B_stop1=fft(B_stop1,n); %对滤去800Hz噪声后的采样信号进行傅立叶变换
FY1_B_stop1=abs( FY_B_stop1);
FY_B_stop2=fft(B_stop2,n); %对滤去800Hz和200Hz噪声后的采样信号进行傅立叶变换
FY1_B_stop2=abs( FY_B_stop2);
FY_B_stop3=fft(B_stop3,n); %对滤去800Hz和200Hz和50Hz噪声后的采样信号进行傅立叶变换
FY1_B_stop3=abs( FY_B_stop3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 画出用不同的滤波器滤去相应噪声频率后的频域比较图 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,1);plot(f, FY1(1:n1+1));grid on;
xlabel('频率（f）');ylabel(' 幅度（FY1）'); title('加噪声后的频谱图');
subplot(2,2,2);plot(f, FY1_B_stop1(1:n1+1));grid on;
xlabel('频率（f）');ylabel(' 幅度（ FY1_B_stop1）'); title('用带阻滤波器滤去800Hz噪声后的的频谱图');
subplot(2,2,3);plot(f, FY1_B_stop2(1:n1+1));grid on;
xlabel('频率（f）');ylabel(' 幅度（ FY1_B_stop2）'); title('用带阻滤波器滤去800Hz和200Hz噪声后的的频谱图');
subplot(2,2,4);plot(f, FY1_B_stop3(1:n1+1));grid on;
xlabel('频率（f）');ylabel(' 幅度（ FY1_B_stop3）'); title('用高通滤波器滤去800Hz和200Hz和50Hz噪声后的的频谱图');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%读出和写出滤波后的语音信号 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pause(5); % 暂停5s
wavplay(B_stop3,fs); % 读出滤去800Hz和200Hz和50Hz噪声后的采样信号
wavwrite(B_stop3,fs,16,'0601331007巴特沃思滤波器滤波后的语音信号.wav');% 写出滤波后的语音信号

切比雪夫I型滤波器滤波

%%%%%%%%%%%%%%%%%%%%%%%%%%% 参数初始化 %%%%%%%%%%%%%%%%%%%%%%

clear all; % 关闭以往仿真的数据栈
A=5; %噪声的振幅
u1=0.01;%噪声振幅的衰减系数
u2=0.02;
u3=0.03;

f1=50;%多频噪声频率
f2=200;
f3=800;

[Y,fs,bits]=wavread('0601331007.wav'); %读出信号1，采样率和采样位数
% sound(Y,fs);
% [X,Fs,bits]=wavread('zhanglow.wav'); %读出信号2，采样率和采样位数
n=length(Y); %求采样信号1的长度
% m=length(X); %求采样信号2的长度
t=0:1/fs:(n-1)/fs; %求采样信号1时域上的采样点数
% t1=0:1/Fs:(m-1)/fs; %求采样信号2时域上的采样点数
% % % % % % t=(0:n-1)/fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 构造噪声 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

y1=u1*A*sin(2*pi*f1*t)+u1*A*sin(2*pi*f2*t)%两两噪声相加
y2=u1*A*sin(2*pi*f3*t)
y=y1+y2
% % % % % % Z1=u1*A*sin(2*pi*f1*t1)+u1*A*sin(2*pi*f2*t1)%两两噪声相加
% % % % % % Z2=u1*A*sin(2*pi*f3*t1)
% % % % % % Z=Z1+Z2

% figure(1) %画出噪声的时域图，取200个点
% % % % % % % % % % plot(t ,y);grid on;
% plot(t(1:200),y(1:200));grid on;
%
% xlabel('时间 (t)');
% ylabel('幅度（y)');
% title('噪声的时域图');
%
% Fy=fft(y,n); %对噪声进行傅立叶变换
% Fy1=abs(Fy);
% n1=floor(n/2); %对zhang求采样点数的一半
% n2=floor(m/2); %对zhanglow求采样点数的一半
% f=(0:n1)*fs/n; %时域上的采样点数
% f1=(0:n2)*Fs/m; %时域上的采样点数
%
% figure(2) %画出噪声的频谱图
% plot(f,Fy1(1:n1+1));grid on;
% xlabel('频率（f）');
% ylabel(' 幅度（Fy1）');
% title('噪声的频谱图');
%
% %%%%%%%%%%%%%%%%%%% 加噪声前后的时域图比较 %%%%%%%%%%%%%%%%%%%%
Y=Y(:,1);%此处假设声音是双声道的，只取单声道作分析
% X=X(:,1);%此处假设声音是双声道的，只取单声道作分析
Y1=y+Y'; %对采样信号进行转置再加上噪声
sound(Y1,fs);
% figure(3) %画出加噪声前后的时域比较图
% subplot(2,1,1);plot(t,Y);grid on;
% plot(t(1:1000),Y(1:1000));grid on;
% xlabel('时间 (t)');
% ylabel('幅度（Y)');
% title('加噪声前的时域图');
%
% subplot(2,1,2);plot(t,Y1);grid on;
% % % % % % % % % % plot(t(1:1000),Y1(1:1000));grid on;
% xlabel('时间 (t)');
% ylabel('幅度（Y1)');
% title('加噪声后的时域图');
%
% %%%%%%%%%%%%%%%%%%%%%%%%%% 加噪声前后的频域图比较 %%%%%%%%%%%%
% FY=fft(Y1,n); %对加噪声后的采样信号1进行傅立叶变换
% FY1=abs(FY);
% FY2=fft(Y,n); %对加噪声前的采样信号1进行傅立叶变换
% FY3=abs(FY2);
% FX=fft(X,m); %对加噪声前的采样信号2进行傅立叶变换
% FX4=abs(FX);
%
% figure(4) %画出加噪声前后的频域比较图
% subplot(2,1,1);
% plot(f,FY3(1:n1+1));grid on;
% xlabel('频率（f）');
% ylabel(' 幅度（FY3）');
% title('加噪声前的频谱图');
%
% subplot(2,1,2);
% plot(f,FY1(1:n1+1));grid on;
% xlabel('频率（f）');
% ylabel(' 幅度（FY1）');
% title('加噪声后的频谱图');
%
% %%%%%%%%%%%%%%%%%%%%%%%% 不同音调相同音量信号的频域、时域图比较 %%%%%%
% figure(5)
% subplot(2,2,1);plot(f,FY3(1:n1+1));grid on;
% xlabel('频率（f）');ylabel(' 幅度（FY3）');title('zhang加噪声前的频谱图');
% subplot(2,2,3);plot(f1,FX4(1:n2+1));grid on;
% xlabel('频率（f）');ylabel(' 幅度（FX4）');title('zhanglow加噪声前的频谱图');
% subplot(2,2,2);plot(t,Y);grid on;
% xlabel('时间 (t)');ylabel('幅度（Y)');title('zhang加噪声前的时域图');
% subplot(2,2,4);plot(t1,X);grid on;
% xlabel('时间 (t1)');ylabel('幅度（X)');title('zhanglow加噪声前的时域图');
% % sound(Y1,fs);%读出加噪后的信号

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 构造切比雪夫I型滤波器 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(6); %画出zhang加噪声后的时域图

plot(t,Y1);
xlabel('时间 (t)');ylabel('幅度（Y1)');title('0601331007加噪声后的时域图');grid on;
% axis([0.04 0.06 -3 3]);

%%%%%%%%%%%%%% 构造800Hz的带阻滤波器 %%%%%%%%%%%%%

f0_stop1=800; %
fc=20; % 设定过渡带的频率
Rp=1; % 通带区的波纹系数
Rs=30; % 阻带区的波纹系数
ws_stop1=[(f0_stop1-fc)/(fs/2) (f0_stop1+fc)/(fs/2)]; % 通带的拐角频率
wp_stop1=[(f0_stop1-5*fc)/(fs/2) (f0_stop1+5*fc)/(fs/2)]; % 阻带的拐角频率
[N_stop1,wn_stop1]=cheb1ord(wp_stop1,ws_stop1,Rp,Rs); % 求出切比雪夫I型滤波器的阶数N及频率参数wc
[num_stop1,den_stop1]=cheby1(N_stop1,Rp,wn_stop1,'stop'); % 求出切比雪夫I型带阻数字滤波器的传递函数模型系数
[h_stop1,w_stop1]=freqz(num_stop1,den_stop1,fs) % 求出离散系统频率响应的数值
figure(7);
subplot(2,3,1); % 画出带阻滤波器的幅频特性图
plot(w_stop1*fs/(2*pi),20*log10(abs(h_stop1)));
xlabel('频率（w_stop1）');ylabel('幅度（h_stop1)');title('带阻滤波器在800Hz处的幅频特性图');grid on;
B_stop1=filter(num_stop1,den_stop1,Y1); %对含噪信号Y1进行带阻滤波

% pause(10);
% sound(B_stop1,fs); % 读出滤去800Hz噪声后的采样信号
subplot(2,3,4); % 画出带阻滤波器滤去800Hz噪声后的时域图
plot(t,B_stop1);
xlabel('时间 (t)）');ylabel('幅度（B_stop1)');title('用带阻滤波器滤去800Hz噪声后的时域图');grid on;
% axis([0.04 0.06 -30 3]);

% %%%%%%%%%%%%%% 构造200Hz的带阻滤波器 %%%%%%%%%%%%%%%%%%%%%%
f0_stop2=200;
ws_stop2=[(f0_stop2-fc)/(fs/2) (f0_stop2+fc)/(fs/2)]; % 通带的拐角频率
wp_stop2=[(f0_stop2-5*fc)/(fs/2) (f0_stop2+5*fc)/(fs/2)]; % 阻带的拐角频率
[N_stop2,wn_stop2]=cheb1ord(wp_stop2,ws_stop2,Rp,Rs); % 求出切比雪夫I型滤波器的阶数N及频率参数wc
[num_stop2,den_stop2]=cheby1(N_stop2,Rp,wn_stop2,'stop'); % 求出切比雪夫I型带阻数字滤波器的传递函数模型系数
[h_stop2,w_stop2]=freqz(num_stop2,den_stop2,fs) % 求出离散系统频率响应的数值

subplot(2,3,2);

plot(w_stop2*fs/(2*pi),20*log10(abs(h_stop2)));
xlabel('频率（w_stop2）');ylabel('幅度（h_stop2)');title('带阻滤波器在200Hz处的幅频特性图');grid on;
B_stop2=filter(num_stop2,den_stop2,B_stop1);%对含噪信号Y1进行带阻滤波

% pause(10);
% sound(B_stop2,fs); % 读出滤去800Hz和200Hz噪声后的采样信号

subplot(2,3,5); % 画出带阻滤波器滤去800Hz和200Hz噪声后的时域图
plot(t,B_stop2);
xlabel('时间 (t)）');ylabel('幅度（B_stop2)');title('带阻滤波器滤去800Hz和200Hz噪声后的时域图');grid on;
% axis([0.04 0.06 -3 3]);
%
%
%%%%%%%%%%%%%% 构造50Hz的带阻滤波器 %%%%%%%%%%%%%%%%%%%%%%
f0_stop3=50
ws_stop3=(f0_stop3-fc)/(fs/2) % 阻带的拐角频率
wp_stop3=(f0_stop3+5*fc)/(fs/2) % 通带的拐角频率
Wws=0.0075;
Wwp=0.0375;

[N_stop3,wc_stop3]=cheb1ord(Wwp,Wws,Rp,Rs,'s'); % 求出切比雪夫I型滤波器的阶数N及频率参数wc
[num_stop3,dem_stop3]=cheby1(N_stop3,Rp,wc_stop3,'high'); % 求出切比雪夫I型滤波器高通数字滤波器的传递函数模型系数
[h_stop3,w_stop3]=freqz(num_stop3,dem_stop3,fs); % 求出离散系统频率响应的数值
subplot(2,3,3); %画出高通滤波器的幅频特性图
plot(w_stop3*fs/(2*pi),(abs(h_stop3)));
xlabel('频率（w_stop3）');ylabel('幅度（h_stop3)');title('高通滤波器在50Hz处的幅频特性图');grid on;
B_stop3=filter(num_stop3,dem_stop3,B_stop2);

subplot(2,3,6); % 画出高通滤波器滤去800Hz和200Hz和50Hz噪声后的时域图
plot(t,B_stop3);
xlabel('时间 (t)）');ylabel('幅度（B_stop3)');title('滤去800Hz和200Hz和50Hz噪声后的时域图');grid on;
% axis([0 10 -60000 60000]);
%%%%%%%%%%%%%%%%%用不同的滤波器滤去相应噪声频率后的时域比较图 %%%%%%%%%%
figure(8);
subplot(2,2,1);plot(t,Y1);grid on;
xlabel('时间 (t)');ylabel('幅度（Y1)');title('加噪声后的时域图');
subplot(2,2,2);plot(t,B_stop1);
xlabel(' 时间 (t)）');ylabel('幅度（B_stop1)');
title('用带阻滤波器滤去800Hz噪声后的时域图');grid on;
subplot(2,2,3);plot(t,B_stop2);
xlabel(' 时间 (t)）');ylabel('幅度（B_stop2)');
title( ' 带阻滤波器滤去800Hz和200Hz噪声后的时域图');grid on;
subplot(2,2,4);plot(t,B_stop3);
xlabel(' 时间 (t)）');ylabel('幅度（B_stop3)');
title( '高通滤波器滤去800Hz