%% 对一段音频滤去噪声
%% Read Portion of Audio File
% Create a FLAC file from the example file |handel.mat|, and then read only
% the first 2 seconds.
%%
% Create a FLAC (|.flac|) file in the current folder.
load handel.mat
filename = 'handel.flac';
audiowrite(filename,y,Fs);
%%
% Read only the first 2 seconds.
samples = [1,1*Fs];
clear y Fs
[y,Fs] = audioread(filename,samples);
%%
% Play the samples.
sound(y,Fs);
pause(1.5);
%% 画y时域图 增加噪声得y2 FFT得Y2
t=(0:length(y)-1)/Fs;
t=t';
noise=0.3*sin(2*pi*1800*t); %1800hz噪声
y2=y+noise;
sound(y2,Fs);
pause(1.5);
figure;plot(t,y);
figure;plot(t,y2);
Y=fft(y);
Y2=fft(y2);
f=(0:length(Y2)-1)/length(Y2) * Fs;
figure;plot(f,abs(Y))
figure;plot(f,abs(Y2))
xlim([0,0.5*Fs])
ylim([0,350])
%% 巴特沃斯带阻滤波器
Fs=8000; %采样频率
fp1=1500;fp2=2100;
fs1=1750;fs2=1850;
wp1=fp1/Fs*2*pi; wp2=fp2/Fs*2*pi; %通带截止频率
wp=[wp1,wp2];
ws1=fs1/Fs*2*pi; ws2=fs2/Fs*2*pi; %阻带截止频率
ws=[ws1,ws2];
Rp=1;As=30;
[n,wc]=buttord(wp/pi,ws/pi,Rp,As)
[b,a]=butter(n,wc,'stop') %求数字带阻滤波器系数
[H,w]=freqz(b,a);
dbH=20*log10(abs(H)/max(abs(H)));
figure;plot(w/2/pi*Fs,dbH,'r');
axis([0,4000,-40,2]);
%% 滤波得y3
y3=filter(b,a,y2);
Y3=fft(y3);
figure;plot(f,abs(Y3));
xlim([0,0.5*Fs])
sound(y3,Fs);