Hello! I have been trying to remove 50 Hz interference and reduce general noise on this ECG result for ages, and just can't get my head around any of the solutions online. Should I use a lowpass filter? How do I implement it? Thanks!
MATLAB: 50 Hz interference and noise reduction from ECG.
interferencenoise
Related Solutions
hello
see below , we are generating a signal with 2 frequencies , then do fft for showing their frequencies.
the time plot will not be very useful but I let it for your info
you can easily adapt it to your own needs
the notch filter is a bonus , to show how filter out one disturbing frequency
hope it helps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% dummy data
Fs = 1000;samples = 25000;t = (0:samples-1)'*1/Fs;signal = cos(2*pi*50*t)+cos(2*pi*100*t)+1e-3*rand(samples,1); % two sine + some noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%NFFT = 1000; %
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');gridtitle(['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(locs+.02,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)'; endfft_spectrum = fft_spectrum(select);freq_vector = (select - 1)*Fs/nfft;end
hello
this code will do fft / spectrogram and bandpass filtering - adapt the corner frequencies to your own needs !
clcclear all%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%NFFT = 4096; %
OVERLAP = 0.5;% 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 = 1;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% data
% [data,Fs] = audioread('myWAVaudiofile.wav');
[data,Fs] = audioread('Immigran_Hiss.wav');channel = 1;signal = data(:,channel);samples = length(signal);%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;if decim>1 signal = decimate(signal,decim); Fs = Fs/decim;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : 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))');endfigure(1),plot(freq,sensor_spectrum_dB,'b');gridtitle(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);xlabel('Frequency (Hz)');ylabel(my_ylabel);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[sg,fsg,tsg] = specgram(signal,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 neededif 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);imagesc(tsg,fsg,sg_dBpeak);colormap('jet');axis('xy');colorbar('vert');gridtitle([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);xlabel('Time (s)');ylabel('Frequency (Hz)');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% removal of hiss between t = 0 and t = 1.68 s
% use of bandpass filter to remove noise below 500 and above 2500 Hz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%data_filtered = data;t_min = 0;t_max = 1.68;ind = 1+fix(t_min*Fs:t_max*Fs-1);% keep only signal content extracted by bandpass filters
f_low = 500;f_high = 2500;N = 4;% signal_filtered = zeros(size(signal));
[b,a] = butter(N,2/Fs*[f_low f_high]);% now let's filter the signal
data_filtered(ind,:) = filtfilt(b,a,data(ind,:));% save as wav file
audiowrite('Immigran_Hiss-filt.wav',data_filtered,Fs);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% display 2B : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[sg,fsg,tsg] = specgram(data_filtered(:,1),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 neededif 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 spectrogramfigure(3);imagesc(tsg,fsg,sg_dBpeak);colormap('jet');axis('xy');colorbar('vert');gridtitle([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);xlabel('Time (s)');ylabel('Frequency (Hz)');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(:));endfunction [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)
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)'; endfft_spectrum = fft_spectrum(select);freq_vector = (select - 1)*Fs/nfft;end
Best Answer