matlab滤音频 某一频率的噪声 利用带阻滤波器

%% 对一段音频滤去噪声


%% 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);