if you look at the noise spectrum, you would see it convers the exact same frequency range as your audio signal
applying a simple low pass filter can remove some of the noise (but certainly not all the noise) , but will also remove some of the original signal content , thus reducing the output fidelity.
the denoising techniques are a bit more complex , either using adaptive filters (if you have a sperate access to the noise source) or other means (a simple serach on internet will show you some tools)
just fyi , you can plot spectrum and spectrogram of your data , and apply filters (adapt to you needs) below :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select audio file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% audiofile = 'original.wav';
audiofile = 'noisy.wav';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select filters to apply
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
option_notch = 0; % 0 = without notch filter , 1 = with notch filter
fc_notch = 4725*[1/3 2/3 1]; % notch center frequencies
p = 0.95; % bandwith parameter (closer to 1 reduces the BW)
option_BPF = 0; % 0 = without bandpass filter , 1 = with bandpass filter
f_low = 20;
f_high = 5000;
N_bpf = 2; % filter order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;
% spectrogram dB scale
spectrogram_dB_scale = 80; % the lowest displayed level is XX dB below the max level
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[signal,Fs] = audioread(audiofile);
[samples,channels] = size(signal);
dt = 1/Fs;
% % select channel (if needed)
% channels = 1:2;
% signal = signal(:,channels);
% [samples,channels] = size(signal);
% time vector
time = (0:samples-1)*dt;
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 0;
if decim>1
Fs = Fs/decim;
for ck = 1:channels
newsignal(:,ck) = decimate(signal(:,ck),decim);
end
signal = newsignal;
end
signal_filtered = signal;
samples = length(signal);
time = (0:samples-1)*1/Fs;
%% 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)]
if option_notch ~= 0
for ck =1:numel(fc_notch)
w0 = 2*pi*fc_notch(ck)/Fs;
% digital notch (IIR)
num1z=[1 -2*cos(w0) 1];
den1z=[1 -2*p*cos(w0) p^2];
% now let's filter the signal
signal_filtered = filter(num1z,den1z,signal_filtered);
end
end
%% band pass filter section %%%%%%
if option_BPF ~= 0
[b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
signal_filtered = filter(b,a,signal_filtered);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1000; %
OVERLAP = 0.75; % percentage of overlap
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
figure(ck),
subplot(211),plot(time,signal(:,ck),'b');grid on
title(['Time plot / Fs = ' num2str(Fs) ' Hz / raw data / Channel : ' num2str(ck)]);
xlabel('Time (s)');ylabel('Amplitude');
subplot(212),plot(time,signal_filtered(:,ck),'r');grid on
title(['Time plot / Fs = ' num2str(Fs) ' Hz / filtered data / Channel : ' num2str(ck)]);
xlabel('Time (s)');ylabel('Amplitude');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);
df = freq(2)-freq(1); % frequency resolution
for ck = 1:channels
figure(channels+ck),
plot(freq,20*log10(spectrum_raw(:,ck)),'b',freq,20*log10(spectrum_filtered(:,ck)),'r');grid on
title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB (L))');
legend('raw','filtered');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
[sg,fsg,tsg] = specgram(signal(:,ck),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));
[sgf,fsgf,tsgf] = specgram(signal_filtered(:,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
sgf_dBpeak = 20*log10(abs(sgf))+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)));
sgf_dBpeak = sgf_dBpeak+(pondA_dB*ones(1,size(sgf_dBpeak,2)));
my_title = ('Spectrogram (dB (A))');
else
my_title = ('Spectrogram (dB (L))');
end
% saturation of the dB range : 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;
sgf_dBpeak(sgf_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(2*channels+ck);
subplot(2,1,1),imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid on
df = fsg(2)-fsg(1); % freq resolution
title([my_title ' Raw Signal / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
xlabel('Time (s)');ylabel('Frequency (Hz)');
subplot(2,1,2),imagesc(tsgf,fsgf,sgf_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid on
title([my_title ' / Filtered Signal / Channel : ' num2str(ck)]);
xlabel('Time (s)');ylabel('Frequency (Hz)');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
