High Frequency Noise baseline filter

I have been dropped head first into a signals project and I have been doing my best to keep my basics about me but I am caught up on a few concepts that some clarification woud be great on.

I have a signal sampled at 80000 Hz and I have the background noise of that signal sampled for 3000 samples before the data trend began ( it is an explosives test so the trigger captured data ~ 5ms before detection).

The background noise is very consistent pre blast. I thought it would be a good idea to remove that noise alone with a bandpass filter first before trying to remove any noise from the experiment itself. But that did not seem to work as I expected as when I took the FFT the Y axis did not make it clear what frequencies to remove (other than the 60Hz for wall noise).

Q1: Is this apporach flawed? It makes sense to me but I hope its just my failure of coding that I can problem solve.

Q2: If it is flawed, what is the best way to remove the underlying sensor noise?

[Jonas]

i would say removing the background noise using the noise profile is a good idea. It would be easier to help if you could provide a sample file and tell us, how exactly you tried to filter (code)?

for a faster preview of what could be the result, you could also at first try to use audactiy, which can do exactly what you were asking for. load the audio into audacity, select the noise, then Effect->Noise Reduction and get noise profile. Then select the rest of the audio and again Effect->Noise Reduction and now you can change the parameters a bit and listen to the resulting sound

[Mathieu NOE]

Problem solved ? otherwise (again) pls share some data

[Kyle]

Attached is one channel of brain data. Sorry for the delay.

[Mathieu NOE]

as you will see in my answer, the noise is presnt over the entire record , not only before the blast.

of course , it's pretty easy to remove the noise by zeroing the signal before the blast, but I wonder what is brings ?

and once the blast occured , the signal amplitude becomes sufficiently large so the background noise shoudn't be a problem anymore.

just my 2 cents

Answer 1

[Star Strider]

Your data have broadband noise, so as might be expected, a frequency-selective filter is going to have very little effect. Probably the most robust approach woudl be to use wavelet denoising (requiring the Wavelet Toolbox) however you can also approach this by using the Savitzky-Golay filter (the sgolayfilt function). I experimented with that here —

load('braindata-1.mat')
% whos
s = double(ParietalDataArrayNOFilt1);
Fs = 80000;
L = numel(s);
t = linspace(0, L-1, L).'/Fs;

figure
plot(t, s)
grid
xlabel('Time')
ylabel('Amplitude')
title('Original')

[FTs1,Fv] = FFT1(s,t);

figure
plot(Fv, abs(FTs1)*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
xlim([0 1E+3])
title('Original')

s_sgf = sgolayfilt(s, 3, 2001);

figure
plot(t, s_sgf)
grid
xlabel('Time')
ylabel('Amplitude')
title('Savitzky-Golay Filtered')

[FTs2,Fv] = FFT1(s_sgf,t);

figure
plot(Fv, abs(FTs2)*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
xlim([0 1E+3])
title('Savitzky-Golay Filtered')

function [FTs1,Fv] = FFT1(s,t)
    s = s(:);
    t = t(:);
    L = numel(t);
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
    Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv);
end

Experiment with the filter arguments to get the result you want.

The Savitzky-Golay filter works essentially as a multiband FIR notch filter, although it works on a different principle. It is extremely useful for signals with broadband noise, since it reduces noise over a relatively wide range of frequencies. I’m not certain what your eventual goal is, however this approach will likely get you closer to a signal with signfiicantly less noise. (One way to eliminate broadband noise as well as mains frequency interference is to use one or more reference sensors to pick up environmental noise and other signals, and then subtract it from the signal sensor using a differential amplifier, ahead of the ADC process. This approach is routinely used in biomedical instrumentation.)

.

[Mathieu NOE]

this is quite a brutal signal filtering

yes, sure , using a low frequency filter will remove a big portion of the noise but from the signal too

the signal peak amplitude is reduced by a factor 10 (also during the useful blast section)

a spectrogram before / after SG filtering shows clearly how the signal has lost a good portion of it's content - is it really what we wanted ?? (why having limited the fft range to 1 kHz as there is significant energy up to 5/10 kHz ?)

I don't see here what the SG filter brings as benefits vs other type of filters , but why not...

anyway, yes, the best option is to have a separate sensor to pick the environmental noise for differential acquisition

Answer 2

[Mathieu NOE]

hello again

an explosive blast record will have a very wide frequency range , the spectrogram will tell you that the signal energy content goes up to 5 to 10 kHz (if I take care of 80 dB of dynamics, counting from the max amplitude to 80 dB lower (or a linear factor of 1 to 1/10,000 in range for max and min values displayed)

the amplitude in dB is given by the color

now the sensor / acquisition system has it's noise spread with constant energy accross the same entire frequency range , so it's not easy to define a fixed parameter filter that can remove most of the noise without significantly distording (reducing) the peak response

the question is how much you can accept of signal peak reduction vs less noise ?

here below some general code I am using to make noise and vibration signals analysis

in the first section you can define optionnal filters to use (typically Butterworth filters , notch filters, you can also use smoothdata which will act as a low pass filter)

NB : I downsampled the data by factor 4, as signal energy above 10 kHz is zero

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select filters to apply (optionnal)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

option_BPF = 1; % 0 = without bandpass filter , 1 = with bandpass filter
    f_low = 20;
    f_high = 5000;
    N_bpf = 3;           % filter order
    
option_smoothdata = 0; % "low pass" filter 

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)

     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('braindata-1.mat')

signal = double(ParietalDataArrayNOFilt1);
clear ParietalDataArrayNOFilt1;

[samples,channels] = size(signal);
Fs = 80e3; % Hz
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 = 4;
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

%% "smoothdata" low pass  filter section %%%%%%
if option_smoothdata ~= 0
    signal_filtered = smoothdata(signal,'gaussian',25);
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