how to remove/reduce power line noise from a short signal?

hello all,

I am trying to reduce / remove power line noise and its harmonic from a short signal. The frequency resolution of the fft of the signal is 40 Hz which makes using the notch filter very challenging, Is there any method that I can apply to reduce powerline noise without hurting the signal iteslef?

[Star Strider]

Powerline noise (typically 50 - 60 Hz) will be aliased into your signal with a 40 Hz sampling frequency (20 Hz Nyquist frequency). It is highly unlikely that you can remove it. You would either have to use a hardware prefilter to remove it from the original signal (a Bessel lowpass anti-aliasing filter with a cutoff (passband) frequency of 20 Hz would work), or increase the sampling frequency to at least 150 Hz in a new recording in order to be able to remove it from the recorded signal.

[Roshanak Yazdanpour]

Thanks for replying, the sampling rate I have is 19200 Hz but the data is cut in epochs of 480 samples(non continous) and as such makes the frequency bin size 40 Hz.

[Mathieu NOE]

hello

must be first clarified if 40 Hz is actually the sampling rate, which is not what is described in first instance

a fft resolution of 40 Hz can be the result of a sufficiently high sampling rate, but then used with a lower nfft buffer length.

but why are we speaking fft here as the issue is about doing time data notch filtering ?

[Roshanak Yazdanpour]

So if I use a notch filter , I have to kill a wide range of frequencies, right?

[Star Strider]

No. See my Answer.

First, filter your signal.

Second, do whatever analyses you want to on it.

[Mathieu NOE]

sorry , no

maybe that demo code will help you

this includes one notch filter to remove a 50 hz tone

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

%% dummy data
Fs = 1000;
samples = 5000;
t = (0:samples-1)'*1/Fs;
signal = cos(2*pi*50*t)+cos(2*pi*440*t)+1e-3*rand(samples,1); % two sine + some noise


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = Fs;    % so delta f = 1 Hz
Overlap = 0.75;
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.


%% notch filter section %%%%%%
% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)

fc = 50;       % notch freq
wc = 2*pi*fc;
Q = 5;         % adjust Q factor for wider (low Q) / narrower (high Q) notch
                % at f = fc the filter has gain = 0

w0 = 2*pi*fc/Fs;
alpha = sin(w0)/(2*Q);

            b0 =   1;
            b1 =  -2*cos(w0);
            b2 =   1;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
            
% analog notch (for info)
num1=[1/wc^2 0 1];
den1=[1/wc^2 1/(wc*Q) 1];

% digital notch (for info)
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];

freq = linspace(fc-1,fc+1,200);
[g1,p1] = bode(num1,den1,2*pi*freq);
g1db = 20*log10(g1);

[gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
gd1db = 20*log10(gd1);

figure(1);
plot(freq,g1db,'b',freq,gd1db,'+r');
title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
legend('analog','digital');
xlabel('Fréquence (Hz)');
ylabel(' dB')
 
% now let's filter the signal 
signal_filtered = filtfilt(num1z,den1z,signal);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq,fft_spectrum] = myfft_peak(signal, Fs, NFFT, Overlap);
sensor_spectrum_dB = 20*log10(fft_spectrum);% convert to dB scale (ref = 1)

% demo findpeaks
df = mean(diff(freq));
[pks,locs]= findpeaks(sensor_spectrum_dB,'SortStr','descend','NPeaks',2);

[freq,fft_spectrum_filtered] = myfft_peak(signal_filtered, Fs, NFFT, Overlap);
sensor_spectrum_filtered_dB = 20*log10(fft_spectrum_filtered);% convert to dB scale (ref = 1)

figure(2),semilogx(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before notch filter','after notch filter');
xlabel('Frequency (Hz)');ylabel(' dB')
text(freq(locs)+1,pks,num2str(freq(locs)))

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 overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    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;
        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

[Roshanak Yazdanpour]

Thanks for the message and code, the signal you created has enough length (5000 points sampled at 1000 Hz). This makes the notch filter easy. The data I am working with has been sampled at 19200 Hz and the length of data I have for filtering is only 480 samples. So inorder to design a notch filter and due to frequency bin (40 Hz) I should remove all frequencies between 40 and 80 Hz which is not desireable.

Answer 1

[Star Strider]

‘... the sampling rate I have is 19200 Hz but the data is cut in epochs of 480 samples(non continous) and as such makes the frequency bin size 40 Hz.’

That is entirely different. You can use the approach in How to apply three notch filters of 3rd order? to create your filter. It will filter the original 50 Hz mains noise and selected harmonics harmonics of it. It is straightforward to adapt it to your problem.

The sampling intervals must be constant for this (or any signal processing) approach to work. If they are not, use the Signal Processing Toolbox resample function to interpolate your signal to constant sampling intervals. Then, design the filter to work with your selected sampling frequency in the resampled signal vector.

[Roshanak Yazdanpour]

thanks for your help , I will look at at the code and update you if I had any luck.

[Star Strider]

My pleasure!

I will need more information (preferably the signal as well) to adapt this filter to your needs.