Filtering unwanted noise in Audio ECG Signal

I am having a problem removing 2 instances of noise in a audio file of an ECG,

My main goal is to try and recreate the clean signal in the organge.

The image below shows the frequency spectrum of the signal i am trying to filter

i have been able to filter the noise in the 30HZ to 40Hz range, my struggle is with the large spike at the 5Hz range, ive only partly been able to filter it out, but output signal is not clear enough and still has noise.

Below ive include the code ive used to get this far,

Any help would be greatly appreciated

[Y_Noisy, FsNoisy] = audioread("ECG_Waves\ecg80.wav");

[Y_Clean, FsClean] = audioread("ECG_Waves\ecgclean.wav");

t_Noisy = (1:4000)/FsNoisy;

plot(t_Noisy,Y_Noisy)

title('Noisy ECG plot')

ylabel('Amplitude')

xlabel('Time(s)')

t_Clean = (1:4000)/FsClean;

plot(t_Clean,Y_Clean)

title('Clean ECG plot')

ylabel('Amplitude')

xlabel('Time(s)')

% Part 1

%{

b = max(Y_Noisy)

c = max(Y_Clean)

d = min(Y_Clean)

e = min(Y_Noisy)

%}

%Highest number Noisy

max_Noisy = Y_Noisy(1);

for a = 2:4000

if Y_Noisy(a) > max_Noisy

max_Noisy = Y_Noisy(a);

end

end

%Lowest number Clean

min_Noisy = Y_Noisy(1);

for a = 2:4000

if Y_Noisy(a) < min_Noisy

min_Noisy = Y_Noisy(a);

end

end

%Highest number Clean

max_Clean = Y_Clean(1);

for a = 2:4000

if Y_Clean(a) > max_Clean

max_Clean = Y_Clean(a);

end

end

%Lowest number Clean

min_Clean = Y_Clean(1);

for a = 2:4000

if Y_Clean(a) < min_Clean

min_Clean = Y_Clean(a);

end

end

% Part 2

Mean_N = mean(Y_Noisy)

Mean_C = mean(Y_Clean)

for e = 1:4000

DC_Removed_Noisy(e) = Y_Noisy(e) - Mean_N;

DC_Removed_Clean(e) = Y_Clean(e) - Mean_C;

end

%Mean_N_DCr = mean(DC_Removed_Noisy)

%Mean_C_DCr = mean(DC_Removed_Clean)

for r = 1:1000

thou_noisy(r) = Y_Noisy(r);

thou_clean(r) = Y_Clean(r);

thou_noisy_DC_Removed(r) = DC_Removed_Noisy(r);

thou_clean_DC_Removed(r) = DC_Removed_Clean(r);

end

figure

t_Noisy_thou = (1:1000)/FsNoisy;

plot(t_Noisy_thou,thou_noisy)

hold on

t_Noisy_thou_dc_removed = (1:1000)/FsNoisy;

plot(t_Noisy_thou_dc_removed,thou_noisy_DC_Removed)

title('Noisy ECG plot with DC removed for Comparison')

ylabel('Amplitude')

xlabel('Time(s)')

figure

t_Clean_thou = (1:1000)/FsClean;

plot(t_Clean_thou,thou_clean)

hold on

t_Clean_thou_dc_removed = (1:1000)/FsClean;

plot(t_Clean_thou_dc_removed,thou_clean_DC_Removed)

title('Clean ECG plot with DC removed for comparison')

ylabel('Amplitude')

xlabel('Time(s)')

% PArt 3

for i = 1:1024

FFT_Noisy_DC_Removed(i) = DC_Removed_Noisy(i);

FFT_Clean_DC_Removed(i) = DC_Removed_Clean(i);

end

S1 = fft(FFT_Noisy_DC_Removed ,1024);

f1 = (0:1023)*((FsNoisy/2)/1024);

figure

plot(f1, abs(S1(1:1024)))

xlabel ('Frequency(Hz)');

the ylabel ('Magnitude of Fourier Transform');

S2 = fft(FFT_Clean_DC_Removed ,1024);

f2 = (0:1023)*((FsClean/2)/1024);

figure

plot(f2, abs(S2(1:1024)))

xlabel ('Frequency(Hz)');

ylabel ('Magnitude of Fourier Transform');

% Part 4

[a,b] = butter(8,0.5,'low');

[H,f] = freqz(a,b,1000);

plot (f, db(abs(H))); % plots the magnitude of the frequency response in dB

xlabel('rads/samples'); % labels x-axis of frequency units

ylabel ('Freq Response Magnitude (dB)'); % labels y-axis

sf = filter(a,b,thou_noisy_DC_Removed); % filters signal s with filter [a, b]

% and places filtered signal in array sf

f4 = (0:999)*((FsNoisy/2)/1000);

SF = fft(sf, 1000);

figure

plot(f4, abs(SF(1:1000))); % plots FFT amplitude spectra of filtered and

% unfiltered signals

xlabel ('Frequency(Hz)'); % labels x-axis of spectral plot

ylabel ('Magnitude of Fourier Transform'); % labels y-axis of spectral plot

%bandstop(sf,[48 52], 1000)

% wave = bandstop(sf,[10 50], 1000);

[a1, b1] = bandstop(sf,[45 60],1000,"ImpulseResponse","iir","Steepness",0.5 );

wave = filter(b1,sf);

figure

plot(t_Clean_thou_dc_removed,wave)

ylabel('Amplitude')

xlabel('Time(s)')

wave_fft = fft(wave, 1000);

figure

plot(f4, abs(wave_fft(1:1000))); % plots FFT amplitude spectra of filtered and

% unfiltered signals

xlabel ('Frequency(Hz)'); % labels x-axis of spectral plot

ylabel ('Magnitude of Fourier Transform'); % labels y-axis of spectral plot

Answer 1

[Star Strider]

The 5 Hz signal is well within the spectrum of the normal EKG (0-100 Hz), so a bandstop or other frequency-selective filter is inappropriate, as it will remove some of the EKG energy. The best option is to use a Savitzky-Golay filter (sgolayfilt) to eliminate most of the noise while keeping the essential parts of the EKG signal.

.