dsp.m

这是本人用Matlab实现音频噪声滤除的试验

%1)time domain analysis 
[y,fs] = wavread('assA04.wav'); % read signal 'assAo4.wav'
r1=length(y);       % length of each samples 
t=(0:r1-1)/fs;        % use all samples multiply each sampling period, and then get time of the signal 
figure(1)      
plot(t,y);  % plot the signal in the time domain t
xlabel('time (s)')
ylabel('Amplitude')
Title('Original signal in the time domain')
sound(y);      % Convert vector into sound

%2)Frequency domain analysis
x= fft(y);     % to do the fft convert
m=abs(x);      % get the absolute value and complex magnitude of y
u=(0:(r1/2-1))/((r1/2)-1)*(fs/2);      % take the half frequency from the frequency domain 
l=20*log10(m);   % change m to dB form
figure(2);         
plot(u,l(1:r1/2))   % plot the m in dB form against u which is half of fs
xlabel('Frequency (Hz)')
ylabel('Magnitude in dB form')
Title('Original signal in the frequency domain')
figure(3)
psd(y,fs)
Title('power spectral density of original signal')
figure(4)
specgram(y,256,fs);
colorbar;
Title('spectrogram of origninal signal')
xlabel('time(s)')
ylabel('Frequency (Hz)')

%3)band pass filter
b = fir1(273, [500/(fs/2),2100/(fs/2)],hanning(273+1));   % design a filter which is 273th order bandwith is 500Hz--2100Hz,with a hamming window function bandpass filter
[h, n] = impz(b, 1);	% Compute impulse response
figure(5);
plot(n, h)	% Plot impulse response
xlabel('n')
ylabel('h(n)')
Title('Filter Impulse Response')
[H, w] = freqz(b, 1);	% Compute filter frequency response
figure(6);
plot(w*fs/(2*pi), abs(H))	% Plot filter frequency response
xlabel('Frequency (Hz)')
ylabel('Magnitude')
Title('Filter frequency Response')
z = filter(b, 1, y);	% Compute output signal
figure(7)
plot(t,z)
xlabel('Time (s)')
ylabel('Amplitude')
Title('Output Signal')
r2=length(z);
Y = fft(z);	% Compute DFT of output signal
f =(0:(r2/2-1))/((r2/2)-1)*(fs/2);
m2=abs(Y);
g = 20*log10(m2);
figure(8);
plot(f,g(1:r2/2));	% Plot output magnitude spectrum
xlabel('Frequency (Hz)')
ylabel('Magnitude')
Title('Output Magnitude Spectrum')
figure(9);
psd(z,fs)
Title('power spectral density of output signal')
figure(10);
specgram(z,256,fs);
colorbar;
Title('spectrogram of origninal signal')
sound(100*z,fs);
%--------------------------------------------------------------------------
% compare with hamming and hanning
b = fir1(273, [500/(fs/2),2100/(fs/2)],hamming(273+1));
B = fir1(273, [500/(fs/2),2100/(fs/2)],hanning(273+1));
[h, n] = impz(b, 1);
[H, n] = impz(B, 1);
z = filter(b, 1, y);	
Z = filter(B, 1, y);
v = fft(z);
V = fft(Z);
m1=abs(v);
m2=abs(V);
output1= 20*log10(m1);
output2= 20*log10(m2);
r2=length(z);
r3=length(Z);
u1=(0:(r2/2-1))/((r2/2)-1)*(fs/2);
u2=(0:(r3/2-1))/((r3/2)-1)*(fs/2); 
figure(11)
subplot(2,1,1)
plot(h)
title('Impulse Response of using hamming window')
subplot(2,1,2)
plot(H)
title('Impulse Response of using hanning window')
figure(12)
freqz(b)
title('Frequency Response of using hamming window')
figure(13)
freqz(B)
title('Frequency Response of using hanning window')
figure(14)
subplot(2,1,1)
plot(z)
title('output signal with hamming window')
subplot(2,1,2)
plot(Z)
title('output signal with hanning window')
xlabel('time(s)')
ylabel('Amplitude')
figure(15)
subplot(2,1,1)
plot(u1,output1(1:r2/2))
title('output signal with hamming window after DFT convertion')
subplot(2,1,2)
plot(u2,output2(1:r3/2))
title('output signal with hanning window after DFT convertion')
figure(16)
subplot(2,1,1)
psd(output1,fs)
title('Power Spectral Density of output signal with hamming window')
subplot(2,1,2)
psd(output2,fs)
title('Power Spectral Density of output signal with hanning window')
figure(17)
subplot(2,1,1)
specgram(z,256,fs);
title('spectrogram of output signal with hamming window')
subplot(2,1,2)
specgram(Z,256,fs);
title('spectrogram of output signal with hanning window')
xlabel('time(s)')
ylabel('Frequency (Hz)')
sound(100*z,fs);pause(5)
sound(100*Z,fs)

% scaling factor
amax=2^31-1;
m=0;
mmax=0;
mwc=0;

[y,fs] = wavread('assA04.wav');
b = fir1(273, [500/(fs/2),2100/(fs/2)],hanning(273+1));
z = filter(b, 1, y);

for m=0:15;
worstcase=sum(abs(b)*(2^15-1)*(2^m));
mwc=fix(worstcase);
  if mwc<=amax
      mmax=m
  end
end
Bn=fix(b*(2^mmax));
figure(18)
freqz(Bn);
xlabel('Frequency (Hz)')
ylabel('Magnitude')
Title('Filter frequency Response with the final scaling factor')
opn =filter(Bn,1,y)./(2^mmax);
figure(19);
subplot(2,1,1)
plot(Z)
subplot(2,1,2)
plot(opn);
sound(100*opn,fs)
wavwrite((100*opn),fs,'assA04_sf.wav');
figure(1)
freqz(b, 1)
xlabel('Frequency (Hz)')
ylabel('Magnitude')
Title('Filter frequency Response')