双音多频信号的滤波
clear all;clc;n=0:204;N=256;k=0:255;FS=8000;
%%%%%%%%%%%%%%合成各个信号%%%%%%%%%%%%
x0=sin(2*pi*941*n/8000)+sin(2*pi*1336*n/8000);
x1=sin(2*pi*697*n/8000)+sin(2*pi*1209*n/8000);
x2=sin(2*pi*697*n/8000)+sin(2*pi*1336*n/8000);
x3=sin(2*pi*697*n/8000)+sin(2*pi*1477*n/8000);
x4=sin(2*pi*770*n/8000)+sin(2*pi*1209*n/8000);
x5=sin(2*pi*770*n/8000)+sin(2*pi*1336*n/8000);
x6=sin(2*pi*770*n/8000)+sin(2*pi*1477*n/8000);
x7=sin(2*pi*852*n/8000)+sin(2*pi*1209*n/8000);
x8=sin(2*pi*852*n/8000)+sin(2*pi*1336*n/8000);
x9=sin(2*pi*852*n/8000)+sin(2*pi*1477*n/8000);
%%%%%%%%%%%%%%%%%画出合成的各个信号的时域图%%%%%%%%%%%%%%%%%%%
subplot(4,3,1);stem(n,x1);xlabel('n');ylabel('x1');axis([0 210 -2 2]);grid on;
subplot(4,3,2);stem(n,x2);xlabel('n');ylabel('x2');axis([0 210 -2 2]);grid on;
subplot(4,3,3);stem(n,x3);xlabel('n');ylabel('x3');axis([0 210 -2 2]);grid on;
subplot(4,3,4);stem(n,x4);xlabel('n');ylabel('x4');axis([0 210 -2 2]);grid on;
subplot(4,3,5);stem(n,x5);xlabel('n');ylabel('x5');axis([0 210 -2 2]);grid on;
subplot(4,3,6);stem(n,x6);xlabel('n');ylabel('x6');axis([0 210 -2 2]);grid on;
subplot(4,3,7);stem(n,x7);xlabel('n');ylabel('x7');axis([0 210 -2 2]);grid on;
subplot(4,3,8);stem(n,x8);xlabel('n');ylabel('x8');axis([0 210 -2 2]);grid on;
subplot(4,3,9);stem(n,x9);xlabel('n');ylabel('x9');axis([0 210 -2 2]);grid on;
subplot(4,3,11);stem(n,x0);xlabel('n');ylabel('x0');axis([0 210 -2 2]);grid on;
%%%%%%%%%%%%%%%%%%求出各信号的傅里叶变换%%%%%%%%%%%%%%%%%%%
X0=fft(x0,N);
X1=fft(x1,N);
X2=fft(x2,N);
X3=fft(x3,N);
X4=fft(x4,N);
X5=fft(x5,N);
X6=fft(x6,N);
X7=fft(x7,N);
X8=fft(x8,N);
X9=fft(x9,N);
%%%%%%%%%%%%%%%%%%%画出各信号的频谱图%%%%%%%%%%%%%%%%%%%%%
figure;
subplot(4,3,1);stem(k,X1);xlabel('K');ylabel('X1');axis([0 300 -100 100]);grid on;
subplot(4,3,2);stem(k,X2);xlabel('K');ylabel('X2');axis([0 300 -100 100]);grid on;
subplot(4,3,3);stem(k,X3);xlabel('K');ylabel('X3');axis([0 300 -100 100]);grid on;
subplot(4,3,4);stem(k,X4);xlabel('K');ylabel('X4');axis([0 300 -100 100]);grid on;
subplot(4,3,5);stem(k,X5);xlabel('K');ylabel('X5');axis([0 300 -100 100]);grid on;
subplot(4,3,6);stem(k,X6);xlabel('K');ylabel('X6');axis([0 300 -100 100]);grid on;
subplot(4,3,7);stem(k,X7);xlabel('K');ylabel('X7');axis([0 300 -100 100]);grid on;
subplot(4,3,8);stem(k,X8);xlabel('K');ylabel('X8');axis([0 300 -100 100]);grid on;
subplot(4,3,9);stem(k,X9);xlabel('K');ylabel('X9');axis([0 300 -100 100]);grid on;
subplot(4,3,11);stem(k,X0);xlabel('K');ylabel('X0');axis([0 300 -100 100]);grid on;
%%%%%%%%%%%%%%试听各个信号%%%%%%%%%%%%%%%%%%%
sound(x0,FS);
sound(x1,FS);
sound(x2,FS);
sound(x3,FS);
sound(x4,FS);
sound(x5,FS);
sound(x6,FS);
sound(x7,FS);
sound(x8,FS);
sound(x9,FS);
%sound(X0,FS);
%sound(X1,FS);
%sound(X2,FS);
%sound(X3,FS);
%sound(X4,FS);
%sound(X5,FS);
%sound(X6,FS);
%sound(X7,FS);
%sound(X8,FS);
%sound(X9,FS);
%%%%%%%%%%%%%%%巴特沃兹低通滤波器的设计%%%%%%%%%%%%%
Wp1=2*pi*1000/8000; %归一化通带截止频率
Ws1=2*pi*1200/8000; %归一化阻带截止频率
Ap1=1; %通带最大衰减
As1=15; %阻带最小衰减
[N1,Wc1]=buttord(Wp1,Ws1,Ap1,As1); %N为滤波器的阶数；Wc为截频
[b1,a1]=butter(N1,Wc1); %b和a为设计出的IIR数字低通滤波器H(s)的分子多项式和分母多项式的系数矩阵
%%%滤波器的幅频响应？？？？
%%%%%%%%%%%%%%%%%%巴特沃兹高通滤波器的设计%%%%%%%%%%%%%%%%%%%%%%%%
Wp2=2*pi*1500/8000;
Ws2=2*pi*1200/8000;
Ap1=3;
As2=14;
[N2,Wc2]=buttord(Wp2,Ws2,Ap2,As2);
[b2,a2]=butter(N2,Wc2,'high');
%滤波器的幅频响应？？？？？
%%%%%%%%%%%%%%%%%%%%%%%各信号通过低通滤波器后的信号%%%%%%%%%%%%%%%%%%
y01=filter(b1,a1,x0);
y11=filter(b1,a1,x1);
y21=filter(b1,a1,x2);
y31=filter(b1,a1,x3);
y41=filter(b1,a1,x4);
y51=filter(b1,a1,x5);
y61=filter(b1,a1,x6);
y71=filter(b1,a1,x7);
y81=filter(b1,a1,x8);
y91=filter(b1,a1,x9);
%%%%%%%%%%%%%%%%%%%%%%%各信号通过高通滤波器后的信号%%%%%%%%%%%%%%%%%%
y02=filter(b2,a2,x0);
y12=filter(b2,a2,x1);
y22=filter(b2,a2,x2);
y32=filter(b2,a2,x3);
y42=filter(b2,a2,x4);
y52=filter(b2,a2,x5);
y62=filter(b2,a2,x6);
y72=filter(b2,a2,x7);
y82=filter(b2,a2,x8);
y92=filter(b2,a2,x9);
%%%%%%%%%%%%信号通过低通滤波器后的频响、幅频响应%%%%%%%%%%%%%%%%%%%%%%%%%%
Y01=fft(y01);Y01_abs=abs(Y01);Y01_angle=angle(Y01);
Y11=fft(y11);Y11_abs=abs(Y11);Y11_angle=angle(Y11);
Y21=fft(y21);Y21_abs=abs(Y21);Y21_angle=angle(Y21);
Y31=fft(y31);Y31_abs=abs(Y31);Y31_angle=angle(Y31);
Y41=fft(y41);Y41_abs=abs(Y41);Y41_angle=angle(Y41);
Y51=fft(y51);Y51_abs=abs(Y51);Y51_angle=angle(Y51);
Y61=fft(y61);Y61_abs=abs(Y61);Y61_angle=angle(Y61);
Y71=fft(y71);Y71_abs=abs(Y71);Y71_angle=angle(Y71);
Y81=fft(y81);Y81_abs=abs(Y81);Y81_angle=angle(Y81);
Y91=fft(y91);Y91_abs=abs(Y91);Y91_angle=angle(Y91);
%%%%%%%%%%%%信号通过高通滤波器后的频响、幅频响应%%%%%%%%%%%%%%%%%%%%%%%%%%
Y02=fft(y02);Y02_abs=abs(Y02);Y02_angle=angle(Y02);
Y12=fft(y12);Y12_abs=abs(Y12);Y12_angle=angle(Y12);
Y22=fft(y22);Y22_abs=abs(Y22);Y22_angle=angle(Y22);
Y32=fft(y32);Y32_abs=abs(Y32);Y32_angle=angle(Y32);
Y42=fft(y42);Y42_abs=abs(Y42);Y42_angle=angle(Y42);
Y52=fft(y52);Y52_abs=abs(Y52);Y52_angle=angle(Y52);
Y62=fft(y62);Y62_abs=abs(Y62);Y62_angle=angle(Y62);
Y72=fft(y72);Y72_abs=abs(Y72);Y72_angle=angle(Y72);
Y82=fft(y82);Y82_abs=abs(Y82);Y82_angle=angle(Y82);
Y92=fft(y92);Y92_abs=abs(Y92);Y92_angle=angle(Y92);
%%%%低通滤波器的幅度响应图%%%
figure;
subplot(4,3,1);stem(k,Y11_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,2);stem(k,Y21_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,3);stem(k,Y31_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,4);stem(k,Y41_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,5);stem(k,Y51_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,6);stem(k,Y61_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,7);stem(k,Y71_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,8);stem(k,Y81_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,9);stem(k,Y91_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,11);stem(k,Y01_abs);xlabel('{\omega}');ylabel('幅度');grid on;
%%%%低通滤波器的频率响应图%%%
figure;
subplot(4,3,1);stem(k,Y11_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,2);stem(k,Y21_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,3);stem(k,Y31_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,4);stem(k,Y41_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,5);stem(k,Y51_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,6);stem(k,Y61_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,7);stem(k,Y71_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,8);stem(k,Y81_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,9);stem(k,Y91_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,11);stem(k,Y01_angle);xlabel('{\omega}');ylabel('相位');grid on;
%%%%高通滤波器的幅度响应图%%%
figure;
subplot(4,3,1);stem(k,Y12_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,2);stem(k,Y22_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,3);stem(k,Y32_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,4);stem(k,Y42_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,5);stem(k,Y52_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,6);stem(k,Y62_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,7);stem(k,Y72_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,8);stem(k,Y82_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,9);stem(k,Y92_abs);xlabel('{\omega}');ylabel('幅度');grid on;
subplot(4,3,11);stem(k,Y02_abs);xlabel('{\omega}');ylabel('幅度');grid on;
%%%%高通滤波器的频率响应图%%%
figure;
subplot(4,3,1);stem(k,Y12_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,2);stem(k,Y22_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,3);stem(k,Y32_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,4);stem(k,Y42_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,5);stem(k,Y52_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,6);stem(k,Y62_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,7);stem(k,Y72_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,8);stem(k,Y82_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,9);stem(k,Y92_angle);xlabel('{\omega}');ylabel('相位');grid on;
subplot(4,3,11);stem(k,Y02_angle);xlabel('{\omega}');ylabel('相位');grid on;
