clc; % 清除工作区
clear; % 清除命令窗口
close all; % 关闭窗口
%% 读取音频文件
[x,fs]=audioread('audio.mp3'); % 读取音频信号：x是数据，fs是采样率
x = x(:,1)'; % 取第一列数据转置
N=length(x); % 计算信号x的长度
t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率
% sound(x,fs); % 播放音频
X=abs(fft(x));X=X(1:N/2); % 对原始信号进行ff变换，截取幅度谱前半部分
deltaf=fs/N; % 计算频谱的谱线间隔
f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围
figure(1)
subplot(2,1,1);plot(t,x);axis([0 6 -1 1]);grid on; % 原始语音信号的时域图
xlabel('时间/s','FontSize',16);ylabel('幅度/v','FontSize',16);
title('原始语音信号','FontSize',20);
subplot(2,1,2);plot(f,X);axis([0 15000 0 1200]);grid on; % 画原始语音信号幅度谱图
xlabel('频率/Hz','FontSize',16);ylabel('|X(f)|','FontSize',16);
title('原始语音信号幅度谱图','FontSize',20);
%% 生成窄带高斯白噪声
noise = wgn(1,N,-20);
[n,Wn,beta,ftype] = kaiserord([10000 10100 10900 11000],...
[0 1 0],[0.01 0.02 0.01],fs); % 估计kaiser窗的参数
n=n+rem(n,2); % 保证滤波器系数长N+1为奇数
b_noise=fir1(n,Wn,ftype,kaiser(n+1,beta)); % kaiser窗函数滤波器设计
noise_fil=fftfilt(b_noise,noise); % 滤出窄带噪声
%% 对原始信号加窄带高斯白噪声
y=x+noise_fil; % 将窄带高斯白噪声叠加到原始音频上
% sound(y,fs); % 播放加噪信号
Y=abs(fft(y));Y=Y(1:N/2); % 对加噪信号进行ff变换，截取幅度谱前半部分
figure(2)
subplot(2,1,1);plot(t,y);axis([0 6 -1 1]);grid on; % 加噪语音信号的时域图
xlabel('时间/s','FontSize',16);ylabel('幅度/v','FontSize',16);
title('加入干扰后的语音信号','FontSize',20);
subplot(2,1,2);plot(f,Y);axis([0 15000 0 1200]);grid on; % 加噪语音信号频谱
xlabel('频率/Hz','FontSize',16);ylabel('|Y(f)|','FontSize',16);
title('加入干扰后的语音信号幅度谱图','FontSize',20);
%% 滤波器设计
Fs2=fs/2; % 奈奎斯特频率
fp=9900; % 通带频率
fr=10000; % 阻带频率
Rp=3; % 通带最大波纹衰减
As=40; % 阻带最小衰减
delta1 = (10^(Rp/20)-1)/(10^(Rp/20)+1); % 通带波纹线性值
delta2 = (1+delta1)*(10^(-As/20)); % 阻带衰减线性值
ff=[fp fr]/Fs2; % 频率指标（归一化频率）
A=[1 0]; % 和幅值指标
dev=[delta1 delta2]; % 偏离指标（偏离程度）
[M,Wn,beta,ftype] = kaiserord(ff,A,dev); % 估计kaiser窗的参数
M=M+rem(M,2); % 保证滤波器系数长N+1为奇数
n=0:M; % 时间范围
b=fir1(M,Wn,ftype,kaiser(M+1,beta)); % kaiser窗函数滤波器设计
[db,mag,pha,grd,w]=freqz_m(b,1); % 计算频率响应、相位响应、脉冲响应
figure(3)
subplot(2,2,1);plot(w/pi,db);axis([0 1 -150 10]);grid on; % 幅度响应(DB)
xlabel('w/pi','FontSize',16);ylabel('dB','FontSize',16);
title('滤波器幅度响应图','FontSize',20);
subplot(2,2,2);plot(w/pi,mag);axis([0 1 -0.2 1.2]);grid on; % 幅度响应
xlabel('wp','FontSize',16);ylabel('幅度 mag','FontSize',16);
title('滤波器幅度响应图','FontSize',20);
subplot(2,2,3);plot(w/pi,pha);grid on; % 相位响应
xlabel('w/pi','FontSize',16);ylabel('相位 pha','FontSize',16);
title('滤波器相位响应图','FontSize',20);
subplot(2,2,4);stem(n,b);grid on; % 脉冲响应
xlabel('n','FontSize',16);ylabel('h(n)','FontSize',16);
title('滤波器脉冲响应图','FontSize',20);
%% 信号滤波处理
y_fil=fftfilt(b,y); % 用设计好的滤波器对y进行滤波
Y_fil=abs(fft(y_fil)); Y_fil=Y_fil(1:N/2); % 计算频谱取前一半
figure(4)
subplot(3,2,1);plot(t,x);axis([0 6 -1 1]);grid on;
xlabel('时间/s','FontSize',16);ylabel('幅度/v','FontSize',16);
title('原始语音信号','FontSize',20);
subplot(3,2,2);plot(f,X);axis([0 15000 0 1200]);grid on;
xlabel('频率/Hz','FontSize',16);ylabel('|X|(f)','FontSize',16);
title('语音信号幅度谱','FontSize',20);
subplot(3,2,3);plot(t,y);axis([0 6 -1 1]);grid on;
xlabel('时间/s','FontSize',16);ylabel('幅度/v','FontSize',16);
title('加干扰后的语音信号','FontSize',20);
subplot(3,2,4);plot(f,Y);axis([0 15000 0 1200]);grid on;
xlabel('频率/Hz','FontSize',16);ylabel('|Y|(f)','FontSize',16);
title('加干扰语音信号幅度谱','FontSize',20);
subplot(3,2,5);plot(t,y_fil);axis([0 6 -1 1]);grid on;
xlabel('时间/s','FontSize',16);ylabel('幅度/v','FontSize',16);
title('滤波后语音信号','FontSize',20);
subplot(3,2,6);plot(f,Y_fil);axis([0 15000 0 1200]);grid on;
xlabel('频率/Hz','FontSize',16);ylabel('|Y\_fil|(f)','FontSize',16);
title('滤波后语音信号幅度谱','FontSize',20);
% sound(y_fil,fs); % 播放音频
