MATLAB Answers

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

35 views (last 30 days)
Roshanak Yazdanpour
Roshanak Yazdanpour on 22 Jan 2021
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?

  7 Comments

Show 4 older comments
Star Strider
Star Strider on 22 Jan 2021
No. See my Answer.
First, filter your signal.
Second, do whatever analyses you want to on it.
Mathieu NOE
Mathieu NOE on 22 Jan 2021
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
Roshanak Yazdanpour on 22 Jan 2021
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.

Sign in to comment.

Answers (1)

Star Strider
Star Strider on 22 Jan 2021
... 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.

  2 Comments

Roshanak Yazdanpour
Roshanak Yazdanpour on 22 Jan 2021
thanks for your help , I will look at at the code and update you if I had any luck.
Star Strider
Star Strider on 22 Jan 2021
My pleasure!
I will need more information (preferably the signal as well) to adapt this filter to your needs.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!