% [z,Fs]=audioread('Audio.wav');
% [z,Fs]=audioread('msg.wav');
% [z,Fs]=audioread('shake.wav');
% [z,Fs]=audioread('tweet.wav');
% [z,Fs]=audioread('Global.wav');
[z,Fs]=audioread('system.wav'); %读取语音信号的数据，赋给变量z
z1=rand(length(z),1);                       %产生与z长度一致的均匀分布随机信号

t=(0:(size(z)-1))/Fs; 
A=(min(abs(z))+max(abs(z)))*0.5;
z2=A*sin(2*pi*5000*t);         %加入正弦噪音

z4=z+z1+z2';                    %随机噪声合成
sound(z4);                 %播放滤波前的语音信号

F_t=8000;
F_p=800;
Fs_1=2000;
W_p=2*pi*F_p/F_t;
W_s=2*pi*Fs_1/F_t;
[N_1,Wn_1]=buttord(W_p,W_s,2,17,'s');        %求低通滤波器的阶数和截止频率

[b_1,a_1]=butter(N_1,Wn_1,'s');               %求S域的频率响应的参数
[b_z,a_z]=bilinear(b_1,a_1,0.5); %利用双线性变换实现滤波响应S域到Z域的变换

[H,W]=freqz(b_z,a_z);                %绘制频率响应曲线
figure(1);
plot(W*Fs_1/(2*pi),abs(H));
set(gca,'Fontsize',22);
title('低通滤波器频率响应曲线','Fontsize',22);
xlabel('频率','Fontsize',22);
ylabel('幅值','Fontsize',22);
grid on;

figure(2);
f2=filter(b_z,a_z,z4);
subplot(2,1,1);
plot(t,z4);
set(gca,'Fontsize',22);
title('滤波前加噪声的语音信号时域波形图','Fontsize',22);
xlabel('时间（s）','Fontsize',22);
ylabel('幅值','Fontsize',22);
subplot(2,1,2);
plot(t,f2);                              %画出滤波后的时域图
set(gca,'Fontsize',22);
title('滤波后的语音信号时域波形图','Fontsize',22);
xlabel('时间（s）','Fontsize',22);
ylabel('幅值','Fontsize',22);
sound(f2);                              %播放滤波后的信号
F0=fft(f2,2048);
f=Fs*(0:1023)/2048;
figure(3);
y2=fft(z4,2048);
subplot(2,1,1);
plot(f,abs(y2(1:1024)));                  %画出滤波前的频谱图
set(gca,'Fontsize',22);
title('滤波前加噪声的语音信号的频谱','Fontsize',22)
xlabel('频率','Fontsize',22);
ylabel('幅值','Fontsize',22);
subplot(2,1,2);
plot(f,abs(F0(1:1024)));              %画出滤波后的频谱图
set(gca,'Fontsize',22);
title('滤波后的语音信号的频谱','Fontsize',22);   
xlabel('频率','Fontsize',22);
ylabel('幅值','Fontsize',22);
