音频增加噪声低通滤波降噪的过程

啊

clear all;
close all;
clc;
[s,Fsample,B]=wavread('music.wav');
Ts=1/Fsample;
M=length(s);%获取音频的数据长度
m=0:M-1;

figure;
subplot(3,1,1); plot(m,s)
xlabel('Number of samples'); ylabel('Sample value x(n)');
axis([0 M+100 -1.5 1.5]);grid;


sigma=mean(s.*s);%噪声方差
v= sqrt(sigma)*randn([M,1]);%噪声生成
x= s+v;
wavwrite(x,'music_noise.wav');
f=[0:M/2]*Fsample/M;
%f=Fsample*linspace(0,1,M/2+1);
X=fft(x);
magX=abs(X);


subplot(3,1,2); plot(m,x)
xlabel('Number of samples'); ylabel('Sample value x(n)');
axis([0 M+100 -1.5 1.5]);grid;

subplot(3,1,3);
semilogx(f,magX(1:M/2+1));
xlabel('Frequency (Hz)');
ylabel('Amplitude |X(f)|');
axis([0 max(f)+100*Fsample/M 0 400]);grid;

%窗函数方法计算滤波器
wc=2*pi*2000/Fsample; %截止频率
Nwin = 147;% 查表N= 3.3 *8000/(fs -fp)
hwin=fir1(Nwin-1, wc/pi);%计算低通滤波器系数
[Hwin f1]=freqz(hwin,1,44100,Fsample);%计算低通滤波器的频率响应
magHwin=abs(Hwin);%幅频响应
phaHwin=angle(Hwin);%相频响应
phaHwin=180*unwrap(phaHwin)/pi;
nwin=0:Nwin-1;
figure;
subplot(3,1,1);stem(nwin,hwin);grid;
xlabel('n');ylabel('Impulse response of LPF');
subplot(3,1,2);semilogx(f1,20*log(magHwin));grid;
xlabel('Frequency(Hz)');ylabel('Magnitude response (db)');
subplot(3,1,3);plot(f1,phaHwin);grid;
xlabel('Frequency(Hz)');ylabel('Phase response (degree)');

%滤波之后的信号
y=filter(hwin,1,x);
Y=fft(y);
magY=abs(Y);%滤波之后的频谱

figure;%绘图，滤波后信号的时域和频率波形
subplot(2,1,1);plot(m,y);
xlabel('Number of samples');ylabel('Sample valu y(n)');grid;
subplot(2,1,2);semilogx(f,magY(1:M/2+1));
xlabel('Frequency(Hz)');ylabel('Amplitude |X(f)|');grid;
wavwrite(y,'music_denoise.wav');