Filtering noise from an audio signal

Hi, I've this file and I need to filter it and make it clear I used this code but the voice is still not clear 
Can someine help please?



clc;
close all; 
clear all;

f=0.8;
n=6;
a=fir1(n,f,'high'); % fir high pass filter 
b=fir1(n,f,'low'); % fir law pass filter 
[y,fs] =audioread ('noisy2_group2.wav'); % load asudio file 
%sound(y,fs);
o=filter(a,1,y); % passing audio to designed high pass filter first
p=filter(b,1,o); % passing o to law pass filter 
fvtool(p,1);     % use to display designed fir filter
subplot(2,1,1);   
plot(y);         % Orignal signal y
subplot(2,1,2);
plot(p);         % filter output 

[Walter Roberson]

sound(p, fs)

[Abdullah Adam]

thank you so much but the sound is still not clear, any ideas to make it clear as much as possible?

[Walter Roberson]

Radio astronomers used to have problems with getting slightly fuzzy images. They decided that this meant that they needed to build higher-quality images that were more accurate and less prone to noise, so they did that. The fuzz increased . So they thought they had made a mistake, and built even more carefully, better quality still. The fuzz increased again .

At this point they decided that maybe the fuzz was really there. They studied it, and found that it was what is now known as Cosmic Microwave Background Radiation, the fine details of which turned out to be very important.

So... how do you know that the parts you are trying to get rid of, are not in fact the really important parts of the signal?

[Abdullah Adam]

Ummm how do I know?

[Walter Roberson]

I certainly do not know what the real information is in the signal, so I cannot say anything about how to make it clearer.

I have not done much audio filtering, so I am not the best person to advise you.

[Abdullah Adam]

Thanks for your help

Answer 1

[Mathieu NOE]

hello

my recommendation : define the filters after you have analysed the data (see the spectrogram plot : most of the noise energy lies above 1000 Hz , so I set the cut off frequency of the low pass filter at 1500 hz as a compromise between noise removal vs keep speech intelligible)

here my "optimal" settings ussind a Butterworth high and low pass filters (the notch is not for your case very usefull)

hope it helps

clc
clearvars

option_notch = 0; % 0 = without notch filter , 1 = with notch filter
    fc_notch = 50;       % notch freq

option_LPF = 1; % 0 = without low pass filter  , 1 = with low pass filter
    fc_lpf = 1500;       % LPF cut off freq
    N_lpf = 2;           % filter order

option_HPF = 1; % 0 = without high pass filter  , 1 = with high pass filter
    fc_hpf = 150;       % HPF cut off freq
    N_hpf = 2;           % filter order

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[signal,Fs] = audioread('noisy2_group2.wav');

[samples,channels] = size(signal);

% create time vector
dt = 1/Fs;
time = (0:samples-1)*dt;

%% notch filter section %%%%%%
% y(n)=G*[x(n)-2*cos(w0)*x(n-1)+x(n-2)]+[2*p cos(w0)*y(n-1)-p^2 y(n-2)]
% this difference equation can be converted to IIR filter numerator /
% denominator
if option_notch ~= 0
    w0 = 2*pi*fc_notch/Fs;
    p = 0.995;

    % digital notch (IIR)
    num1z=[1 -2*cos(w0) 1];
    den1z=[1 -2*p*cos(w0) p^2];

    % now let's filter the signal 
    signal = filter(num1z,den1z,signal);
end

%% low pass filter section %%%%%%
if option_LPF ~= 0
    w0_lpf = 2*fc_lpf/Fs;
    % digital notch (IIR)
    [b,a] = butter(N_lpf,w0_lpf);
    % now let's filter the signal 
    signal = filter(b,a,signal);
end

%% high pass filter section %%%%%%
if option_HPF ~= 0
    w0_hpf = 2*fc_hpf/Fs;
    % digital notch (IIR)
    [b,a] = butter(N_hpf,w0_hpf,'high');
    % now let's filter the signal 
    signal = filter(b,a,signal);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 512;    % 
OVERLAP = 0.75;

% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 

% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;

%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),decim);
    Fs = Fs/decim;
    end
   signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;

%%%%%% legend structure %%%%%%%%
for ck = 1:channels
    leg_str{ck} = ['Channel ' num2str(ck) ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),plot(time,signal);grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend(leg_str);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);

% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end


figure(2),plot(freq,sensor_spectrum_dB);grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
legend(leg_str);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
    [sg,fsg,tsg] = specgram(signal(:,ck),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  

    % FFT normalisation and conversion amplitude from linear to dB (peak)
    sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only

    % apply A weigthing if needed
    if option_w == 1
        pondA_dB = pondA_function(fsg);
        sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
        my_title = ('Spectrogram (dB (A))');
    else
        my_title = ('Spectrogram (dB (L))');
    end

    % saturation of the dB range : 
    % saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
    min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
    sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

    % plots spectrogram
    figure(2+ck);
    imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
    axis('xy');colorbar('vert');grid on
    df = fsg(2)-fsg(1); % freq resolution 
    title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Frequency (Hz)');

end

sound(signal,Fs);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end


function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

[Mathieu NOE]

hello

problem solved?

[Abdullah Adam]

Sorry for the late reply, I'll check now