Remove noise on a wav file

Question

I'm working on a small code to learn signal processing on Matlab. I have got a .wav sound with some noise and I just want to remove the noise. I tried the code under but noise is not removed correctly. My idea is to do a cut band filter to remove the different noise components on the fft. After a lot of researches, I don't understant where is my problem. Here my code :

clear all;
clc;
close all;

% We load and plot the signal
[sig,Fs] = audioread('11.wav');

N = length(sig);

figure,plot(sig); title 'signal'

% FFT of the signal
fft_sig = abs(fft(sig)); 

% Normalisation of the frequencies for the plot
k = 0 : length(sig) - 1;  
f = k*Fs/length(sig);

plot(f,fft_sig);

% Loop with 2 elements because there are 2 amplitudes (not sure about
% this)
 for i = 1:2

% I put noise components in an array
[y(i),x(i)] = max(fft_sig);

% Normalisation of the frequencies to eliminate
f(i) =  x(i) * (Fs/N);

% Cut band filter with elimination of f from f-20 to f+20 Hz
b = fir1(2 , 2 * [f(i)-20 f(i)+20] / Fs, 'stop')

sig = filter(b,1,sig);

figure,plot(abs(fft(sig)));title 'fft of the signal'

end

Here the images I got, the fft plot is exactly the same before and after applying the filter, there is a modification only on the x axis :

The sampling frequency is Fs = 22050.

Thanks in advance for your help, I hope I'm clear enough in my description

try using 'freqz' to see what your filter looks like. you want to build a notch filter. if you have the dsp toolbox you can use 'iirnotch' — crowdedComputeeer
It might be an issue with your normalization. fir1 expects units of rad/sec and you have not accounted for this anywhere. You can also discard the second half of the output of fft as it is just a mirror image of the first half (otherwise this may mess up your peak component detection scheme). Regarding your peak component detection, as it stands, your call to max(fft_sig) will return the same exact component every time whereas you need it to output the next largest component on every iteration. — Falimond
And what is the most important - plot the magnitude in logarithmic scale! — jojek

SleuthEye SleuthEye · Accepted Answer · 2015-09-21T17:53:05

Since you haven't explicitly said so, the code you provided essentially defines your noise as narrowband interference at two frequencies (even though that interference may look less 'noisy').

The first thing to notice is that the value x(i) obtained from max(fft_sig) corresponds to the 1-based index of the located maximum. It won't make a huge difference for large N, but it may have an impact for smaller values especially when trying to design a very narrow notch filter. The corresponding frequency would then be:

freq =  (x(i)-1) * (Fs/N);

Also if you are going to iteratively attenuate frequency components you would need to update fft_sig which you use to pick the frequency components to be attenuated (otherwise you will always pick the same component). The simplest would be to recompute fft_sig from the filtered sig:

sig = filter(b,1,sig);
fft_sig = abs(fft(sig));

Finally, since you are trying to attenuate a very narrow frequency range, you may find that you need to increase the FIR filter order by a few orders of magnitudes to get a good attenuation at the desired frequency without attenuating everything else. As was pointed out in comments such a narrow notch filter can often be better implemented using an IIR filter.

The updated code would then look somewhat like:

% Loop with 2 elements because there are 2 amplitudes
for i = 1:2

  % I put noise components in an array
  % Limit peak search to the lower-half of the spectrum since the
  % upper half is symmetric
  [y(i),x(i)] = max(fft_sig(1:floor(N/2)+1));

  % Normalisation of the frequencies to eliminate
  freq =  (x(i)-1) * (Fs/N);

  % Cut band filter with elimination of f from f-20 to f+20 Hz
  b = fir1(2000 , 2 * [freq-20 freq+20] / Fs, 'stop')

  sig = filter(b,1,sig);
  fft_sig = abs(fft(sig));

  %figure;plot(f, abs(fft_sig));title 'fft of the signal'
  figure;plot(f, 20*log10(abs(fft_sig)));title 'fft of the signal'
end

As far as the last FFT plot showing a different x-axis, it is simply because you omitted to provide an x-axis variable (in occurrence f used in the previous plot), so the plot shows the array index (instead of frequencies) as the x-axis.

Remove noise on a wav file

1 Answers