Optimal Cutoff Frequency for Static Noise Detection in ECG Signals?

Hello MATLAB Community,

I am currently working on optimizing the detection of static noise in ECG signals and would greatly appreciate your expertise. Specifically, I am looking to determine the best cutoff frequency for filtering this noise. Below, I have listed the SNR (Signal-to-Noise Ratio) values for different cutoff frequencies in two leads, Lead I and Lead aVL:

• Cutoff Frequency = 0.5 Hz
• SNR in Lead I: 7.98 dB
• SNR in Lead aVL: 5.41 dB
• Cutoff Frequency = 1.0 Hz
• SNR in Lead I: 7.29 dB
• SNR in Lead aVL: 5.11 dB
• Cutoff Frequency = 5.0 Hz
• SNR in Lead I: 4.03 dB
• SNR in Lead aVL: 3.23 dB
• Cutoff Frequency = 10.0 Hz
• SNR in Lead I: 2.17 dB
• SNR in Lead aVL: 1.96 dB

Quantification of Noise:

• Cutoff Frequency = 0.5 Hz
• Number of noise points in Lead I: 299
• Number of noise points in Lead aVL: 341
• Cutoff Frequency = 1.0 Hz
• Number of noise points in Lead I: 278
• Number of noise points in Lead aVL: 304
• Cutoff Frequency = 5.0 Hz
• Number of noise points in Lead I: 179
• Number of noise points in Lead aVL: 213
• Cutoff Frequency = 10.0 Hz
• Number of noise points in Lead I: 127
• Number of noise points in Lead aVL: 137

Additionally, I have attached an image showing the residuals ( lead I ) .

To determine the best threshold value, I used an approach based on minimizing the number of noise points detected in the filtered signal. Here is a detailed explanation of the process:

Process to determine the best threshold value

Calculation of the filtered signal:

For each cutoff frequency, I applied a high-pass filter to remove low-frequency components from the ECG signal.

Threshold definition:

I tested different threshold values, defined as multiples of the standard deviation (STD) of the filtered signal. The tested thresholds were 0.25 * STD, 0.5 * STD, 0.75 * STD, and 1 * STD.

Noise point detection:

For each threshold value, I detected points in the filtered signal where the amplitude exceeds the threshold.

The number of detected noise points is counted for each threshold value.

Selection of the best threshold:

The best threshold is the one that minimizes the number of detected noise points. The hypothesis is that the optimal threshold eliminates noise without affecting the useful components of the ECG signal.

Based on this data, I am seeking advice on the most suitable cutoff frequency for effectively reducing static noise while preserving the integrity of the ECG signal. Any suggestions or insights into methodologies for determining this would be highly valuable.

Thank you in advance for your help!

Best regards,

Answer 1

[Star Strider]

You have not defined ‘static noise’ so I have no idea what you’re actually filtering.

The bandwidth of a normal EKG is generally 0 Hz (D-C) to 45 Hz, although with complex arrythmias (most notably atrial fibrillation), it can go up to 100 Hz, and baseline drift will show up as low-freqauency noise, usually less than 2 Hz. Any lowpass (for only high-frequency noise) or bandpass (for baseline drift with high-frequency noise) filter with cutoff frequencies outside those limits (higher or lower, respectively) will eliminate some information. For isolated mains/powerline frequency noise (usually 50-60 Hz), you can use a notch filter to eliminate it, although a lowpass filter with a 45 Hz cutoff will work as well in that instance.

The filter choice is generally IIR, and I prefer elliptic filters for their computational efficiency.

The best way to design a filter for any purpose is to first calculate the Fourier transform of the signal, see if the noise is band-limited (as opposed to broadband) and then choose the cutoff frequencies with that in mind. For band-limited noise, a frequency-selective filter will work. For broadband noise, I usually use a Savitzky-Golay filter or wavelet denoising, whichever is more efficient for that particular problem.

The signal you posted appears to be relatively ‘clean’, so I would be tempted to just lleave it as it is.

If you want to plot a 2D vectorcardiogram, you can use Lead I and aVF as:

figure

plot(Lead_I, Lead_aVF)

grid

I suggest doing that experiment.

.

[rawaa mejri]

Hello @Star Strider ,

Static noise on ECG data refers to a type of electrical interference that disrupts the recording of the heart's electrical activity. It can manifest as flat lines, sudden spikes, or a crackling pattern superimposed on the ECG signal.

The dataset I am working with (PTB-XL v3) has several issues, including static noise. These issues include signal quality problems such as noise (static_noise and burst_noise), baseline drifts (baseline_drift), and other artifacts such as electrode problems. I need to visualize and then remove these issues, which is why I am currently working on this.

for example :

filename_lr baseline_drift static_noise burst_noise electrodes_problems

records100/21000/21837_lr(empty), I-AVL,(empty)(empty)

Thanks for your help,

[Star Strider]

My pleasure!

I always regarded ‘static noise’ as just noise, as distinguished from mains/powerline noise, that tends to be well-characterised and strictly band-limited.

As I mentioned originally, first calculate the Fourier transform of your signal (also each lead if possible). That will give you the frequency content of the signal as well as the noise, and you can proceed from there to design the filter and determine whatever other signal processing you want or need to do.

[rawaa mejri]

thanks a lot @Star Strider

Here are the Frequency Spectrum for Lead I and aVL , also for all the signal

The best threshold value is : 1.00 * STD

Here is the script ,

% Function to calculate SNR
function snr_value = calculate_snr(signal, noise)
    snr_value = 20 * log10(rms(signal) / rms(noise));
end

folderPath = 'dataset100'; 
fileName = 'records100_21837_lr.csv';
signalFileName = fullfile(folderPath, fileName);

signalData = readtable(signalFileName);

fs = 100; 
t = (0:height(signalData)-1) / fs;

% Extract the signals from the leads
num_leads = size(signalData, 2) - 1; 
all_leads_data = signalData{:, 2:end}; 
combined_signal = sum(all_leads_data, 2); 
combined_fft = fft(combined_signal);
f = (0:length(combined_signal)-1) * fs / length(combined_signal);

figure;
plot(f, abs(combined_fft));
title('Frequency Spectrum of the Combined Signal (All Leads)');
xlabel('Frequency (Hz)');
ylabel('Amplitude');
grid on;

I_data = signalData.I;
AVL_data = signalData.aVL;

figure;
subplot(2,1,1);
I_fft = fft(I_data);
plot(f, abs(I_fft));
title('Frequency Spectrum of Lead I');
xlabel('Frequency (Hz)');
ylabel('Amplitude');
grid on;

subplot(2,1,2);
AVL_fft = fft(AVL_data);
plot(f, abs(AVL_fft));
title('Frequency Spectrum of Lead aVL');
xlabel('Frequency (Hz)');
ylabel('Amplitude');
grid on;

% Test different cutoff frequencies
cutoff_frequencies = [0.5, 1, 5, 10]; % Cutoff frequencies to test

results = struct();

colors = ['r', 'g', 'b', 'm'];

figure;

for i = 1:length(cutoff_frequencies)
    cutoff_freq = cutoff_frequencies(i);
    [b, a] = butter(4, cutoff_freq/(fs/2), 'high');
    
    % Filter the signals to detect static noise
    filtered_I_data = filtfilt(b, a, I_data);
    filtered_AVL_data = filtfilt(b, a, AVL_data);
    
    % Calculate the residues
    noise_I = I_data - filtered_I_data;
    noise_AVL = AVL_data - filtered_AVL_data;
    
    % Calculate the SNR
    snr_I = calculate_snr(I_data, noise_I);
    snr_AVL = calculate_snr(AVL_data, noise_AVL);
    
    results(i).cutoff_freq = cutoff_freq;
    results(i).snr_I = snr_I;
    results(i).snr_AVL = snr_AVL;
    
    subplot(2,1,1);
    hold on;
    plot(t, noise_I, colors(i), 'DisplayName', ['Cutoff = ' num2str(cutoff_freq) ' Hz']);
    hold off;
    
    subplot(2,1,2);
    hold on;
    plot(t, noise_AVL, colors(i), 'DisplayName', ['Cutoff = ' num2str(cutoff_freq) ' Hz']);
    hold off;
end

subplot(2,1,1);
title('Residue Lead I for Different Cutoff Frequencies');
xlabel('Time (s)');
ylabel('Amplitude');
legend show;
grid on;

subplot(2,1,2);
title('Residue Lead aVL for Different Cutoff Frequencies');
xlabel('Time (s)');
ylabel('Amplitude');
legend show;
grid on;

threshold_multipliers = [0.25, 0.5, 0.75, 1]; % Multipliers of the standard deviation

figure;

for i = 1:length(cutoff_frequencies)
    cutoff_freq = cutoff_frequencies(i);
    [b, a] = butter(4, cutoff_freq/(fs/2), 'high');
    filtered_I_data = filtfilt(b, a, I_data);
    filtered_AVL_data = filtfilt(b, a, AVL_data);

    for j = 1:length(threshold_multipliers)
        threshold_I = threshold_multipliers(j) * std(filtered_I_data);
        threshold_AVL = threshold_multipliers(j) * std(filtered_AVL_data);

        static_noise_indices_I = find(abs(filtered_I_data) > threshold_I);
        static_noise_indices_AVL = find(abs(filtered_AVL_data) > threshold_AVL);

        subplot(length(cutoff_frequencies), 2, 2*i-1);
        hold on;
        plot(t, I_data, 'k');
        plot(t(static_noise_indices_I), I_data(static_noise_indices_I), [colors(j) '.'], 'MarkerSize', 10);
        hold off;

        subplot(length(cutoff_frequencies), 2, 2*i);
        hold on;
        plot(t, AVL_data, 'k');
        plot(t(static_noise_indices_AVL), AVL_data(static_noise_indices_AVL), [colors(j) '.'], 'MarkerSize', 10);
        hold off;
    end
end

for i = 1:length(cutoff_frequencies)
    subplot(length(cutoff_frequencies), 2, 2*i-1);
    title(['Lead I Signal with Static Noise for Different Threshold Multipliers (Cutoff = ' num2str(cutoff_frequencies(i)) ' Hz)']);
    xlabel('Time (s)');
    ylabel('Amplitude');
    legend('Original Signal', 'Threshold = 0.25 * STD', 'Threshold = 0.5 * STD', 'Threshold = 0.75 * STD', 'Threshold = 1 * STD', 'Location', 'Best');
    grid on;

    subplot(length(cutoff_frequencies), 2, 2*i);
    title(['Lead aVL Signal with Static Noise for Different Threshold Multipliers (Cutoff = ' num2str(cutoff_frequencies(i)) ' Hz)']);
    xlabel('Time (s)');
    ylabel('Amplitude');
    legend('Original Signal', 'Threshold = 0.25 * STD', 'Threshold = 0.5 * STD', 'Threshold = 0.75 * STD', 'Threshold = 1 * STD', 'Location', 'Best');
    grid on;
end

% Determine the best threshold value
% Use criteria based on SNR and the amount of detected static noise
best_threshold_multiplier = threshold_multipliers(1);
best_num_noise_points = inf;
for i = 1:length(threshold_multipliers)
    threshold_I = threshold_multipliers(i) * std(filtered_I_data);
    static_noise_indices_I = find(abs(filtered_I_data) > threshold_I);
    if length(static_noise_indices_I) < best_num_noise_points
        best_num_noise_points = length(static_noise_indices_I);
        best_threshold_multiplier = threshold_multipliers(i);
    end
end

fprintf('The best threshold value is : %.2f * STD\n', best_threshold_multiplier);

[Star Strider]

Please provide the EKG data as a text file or .mat file (all leads if possible).

[rawaa mejri]

Here is an example(records100_21837.csv).

All the ECG are organised on CSV file ( 13 columns : the first one is the Time , 12 leads ), the attached file contains static noise ( I, aVL leads ).

Thanks a lot.

[Star Strider]

It is important to remove as much noise as possible, while retaining the significant features of the EKG trace for clinical reasons. That is my objective here. Experiment to get the result you want.

My analysis —

T1 = readtable('records100_21837_lr.csv')

T1 = 1000x13 table
    Time      I         II       III       aVR       aVL       aVF        V1        V2        V3        V4        V5        V6  
    ____    ______    ______    ______    ______    ______    ______    ______    ______    ______    ______    ______    ______

       0         0         0         0         0         0         0         0         0         0         0         0         0
    0.01         0    -0.004    -0.004     0.002     0.002    -0.004         0     0.001     0.001     0.008    -0.007    -0.012
    0.02     -0.01    -0.017    -0.007     0.013    -0.001    -0.012         0    -0.006    -0.007         0    -0.021     -0.02
    0.03    -0.017    -0.027    -0.011     0.022    -0.003    -0.019         0    -0.016    -0.022    -0.021     -0.03    -0.011
    0.04    -0.009    -0.023    -0.014     0.016     0.002    -0.019         0    -0.009    -0.016    -0.013    -0.029    -0.009
    0.05     0.007    -0.012    -0.019     0.002     0.013    -0.015     0.001    -0.006    -0.013    -0.009    -0.026    -0.004
    0.06     0.014    -0.003    -0.017    -0.006     0.016     -0.01     0.004    -0.002     -0.01    -0.003    -0.023     0.001
    0.07      0.02     0.002    -0.018    -0.012     0.019    -0.009     0.005     0.004    -0.005     0.004     -0.02    -0.002
    0.08     0.033     0.005    -0.028    -0.019      0.03    -0.012     0.005     0.003    -0.008     0.002    -0.018    -0.006
    0.09     0.035     0.008    -0.026    -0.021      0.03    -0.009     0.005         0     -0.01         0    -0.015    -0.011
     0.1     0.034     0.013    -0.021    -0.024     0.027    -0.004     0.005         0     -0.01         0    -0.014    -0.014
    0.11      0.03     0.015    -0.014    -0.023     0.022         0     0.003         0     -0.01         0    -0.014    -0.017
    0.12     0.029      0.02     -0.01    -0.025      0.02     0.005         0         0     -0.01         0    -0.014    -0.019
    0.13     0.034      0.02    -0.014    -0.027     0.024     0.002    -0.003     0.001     -0.01         0    -0.013    -0.022
    0.14      0.04     0.027    -0.013    -0.034     0.026     0.007    -0.006    -0.012    -0.021    -0.004    -0.019    -0.025
    0.15     0.033     0.017    -0.017    -0.025     0.024         0    -0.009    -0.019    -0.033    -0.009    -0.027     -0.02

VN = T1.Properties.VariableNames;
EKG = T1{:,2:end};

tilenr = [1:2:12 2:2:12];

figure
tiledlayout(6,2)
for k = 1:size(EKG,2)
    nexttile(tilenr(k))
    plot(T1.Time, EKG(:,k))
    grid
    title(VN{k+1})
end
sgtitle('EKG Original')

[FTs1,Fv] = FFT1(EKG,T1.Time);

figure
tiledlayout(6,2)
for k = 1:size(FTs1,2)
    nexttile(tilenr(k))
    plot(Fv, abs(FTs1(:,k)*2))
    grid
    title(VN{k+1})
end
sgtitle('EKG Original — Fourier Transform')

EKGfilt = sgolayfilt(EKG, 4, 7);                           % Savitzky-Golay Filter

figure
tiledlayout(6,2)
for k = 1:size(EKG,2)
    nexttile(tilenr(k))
    plot(T1.Time, EKGfilt(:,k))
    grid
    title(VN{k+1})
end
sgtitle('EKG Filtered')

[FTs2,Fv] = FFT1(EKGfilt,T1.Time);

figure
tiledlayout(6,2)
for k = 1:size(FTs2,2)
    nexttile(tilenr(k))
    plot(Fv, abs(FTs2(:,k)*2))
    grid
    title(VN{k+1})
end
sgtitle('EKG Filtered — Fourier Transform')

figure
tiledlayout(2,2)
nexttile
plot(T1.Time, T1.II)
grid
title('Lead II, Original')
nexttile
plot(T1.Time, EKGfilt(:,2))
grid
ylim([min(ylim) 1])
title('Lead II, Filtered')
nexttile
plot(Fv, abs(FTs1(:,2))*2)
grid
title(["Lead II Original" "Fourier Transform"])
nexttile
plot(Fv, abs(FTs2(:,2))*2)
grid
title(["Lead II Filtered" "Fourier Transform"])
sgtitle('Lead II — Detail')

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

This approach removes some of the high-frequency noise without unduly sacrificing details of the EKG. The problem is that the deflection amplitudes are decreased (some more than others), and since these are important, that is not necessarily a desirable result. Nevertheless, the relative ratios appear to be retained, so multiplying the filtered traces by some appropriate constant would probably be appropriate to restore them to their original amplitudes. I did not do that here, however I would consider it in a clinical setting.

.

[rawaa mejri]

Thank you @Star Strider for your invaluable help as always. My question is, "How can I determine if there are high-frequency noises to eliminate? As I mentioned, the dataset I am working with contains four types of noise (static noise, baseline drift, burst noise, and electrode problems). Regarding 'multiplying the filtered traces by some appropriate constant would probably be appropriate to restore them to their original amplitudes,' I'm not sure how to do that. Do you have any ideas?

Concerning the static noise, is the approach I have taken so far correct? Does it only need filtering, or is it completely wrong?

[Star Strider]

My pleasure.

Use the Fourier transform to see the spectral content of the signal. Then, design the filter appropach appropriatesly.

These signals defy any simple filtering approach, although you could experiment with a multi-band FIR bandstop filter, however the problem with that is that the signals are only 1000 samples in length, and an optimal FIR filter would require a signal length many times that. There is no one good filtering approach with these traces that eliminates most of the noise and retains the essential features so that the filtered EKGs are clinically useful.

Depending on the desired end result, I would experiment with the bandpass function, and simply keep varying the upper stopband (it should be relatively straightforward to fix the lower stopband to eliminate the baseline variation), and then check the resullt for the best feature retention and the least noise. There is no one-best-way to do that.

As for the scaling, use the bounds function to find the minimum and maximul of the unfiltered EKG, and then use the rescale funciton with those values to rescale the filtered EKG. That is the easiest approach.

[rawaa mejri]

Thanks a lot @Star Strider

In the present dataset, there are 4 types of noises (electrodes problems, burst noise, baseline drift, and static noise). Each lead contains some of these noises (e.g., Lead I contains burst noise, baseline drift, or both). Should I remove these noises? Does it affect the signal quality? Is it recommended to leave them in to train the model on data quality similar to reality, resulting in a more robust model?

That would remain unclear to me.

What are your recommendations?

[Star Strider]

As always, my pleasure!

To remove baseline drift, use a highpass filter with an appropriate cutoff frequency (usually about 1 Hz, however thant can vary). Noise with varying frequencies (that I call ‘broadband noise’) can be difficult to deal with. I usually use a Savitzky-Golay filter for that, with varying success, depending principally on the noise amplitude. Band-limited noise can be dealt with using a conventional frequency-based filter (lowpass or bandstop). The most important result is to retain as much of the morpohology (essential features, such as the relative amplitudes of the various deflections, and the essential intervals) as possible, so the result is clinically useful.

I would elimiinate as much noise as possible without significantly changing the essential features of the EKG, and train your routine on whatever is left. Making it robust to various types of noise is important. It may not be possible to elliminate all noise, something that clinical electrographers (as I once was) simply need to deal with.

Important aspects of noise reducetion are correct electrode application, alway using a reference (usually right leg) electrode, always using coaxial cables witth the leads and grouding the shields to the instrument ground, and keeping tthe patient from moving excessively, since that creates baseline drift and muscle activation artifact.

There is only a limited amount you can do other than that.

[rawaa mejri]

Thank you very much @Star Strider.

As I mentioned before, this is the first time I'm working with this kind of data. Could you please help me with leads that have more than one type of noise?

Let's start with the "static noise." How can I detect it (visualize areas with spikes, for example, as shown in the figure in the previous comments)? How do I remove it? Unfortunately, I don't have a specialist with me. Is there a recommended strategy for prioritizing noise removal (for example, in lead I, where there are two different types of noise, how should I remove them)?

[Star Strider]

As always, my pleasure!

It is difficult to isolate the noise in one lead from noise in the other leasds, although that can sometimes occur with bad electrode application.

T1 = readtable('records100_21837_lr.csv')

T1 = 1000x13 table
    Time      I         II       III       aVR       aVL       aVF        V1        V2        V3        V4        V5        V6  
    ____    ______    ______    ______    ______    ______    ______    ______    ______    ______    ______    ______    ______

       0         0         0         0         0         0         0         0         0         0         0         0         0
    0.01         0    -0.004    -0.004     0.002     0.002    -0.004         0     0.001     0.001     0.008    -0.007    -0.012
    0.02     -0.01    -0.017    -0.007     0.013    -0.001    -0.012         0    -0.006    -0.007         0    -0.021     -0.02
    0.03    -0.017    -0.027    -0.011     0.022    -0.003    -0.019         0    -0.016    -0.022    -0.021     -0.03    -0.011
    0.04    -0.009    -0.023    -0.014     0.016     0.002    -0.019         0    -0.009    -0.016    -0.013    -0.029    -0.009
    0.05     0.007    -0.012    -0.019     0.002     0.013    -0.015     0.001    -0.006    -0.013    -0.009    -0.026    -0.004
    0.06     0.014    -0.003    -0.017    -0.006     0.016     -0.01     0.004    -0.002     -0.01    -0.003    -0.023     0.001
    0.07      0.02     0.002    -0.018    -0.012     0.019    -0.009     0.005     0.004    -0.005     0.004     -0.02    -0.002
    0.08     0.033     0.005    -0.028    -0.019      0.03    -0.012     0.005     0.003    -0.008     0.002    -0.018    -0.006
    0.09     0.035     0.008    -0.026    -0.021      0.03    -0.009     0.005         0     -0.01         0    -0.015    -0.011
     0.1     0.034     0.013    -0.021    -0.024     0.027    -0.004     0.005         0     -0.01         0    -0.014    -0.014
    0.11      0.03     0.015    -0.014    -0.023     0.022         0     0.003         0     -0.01         0    -0.014    -0.017
    0.12     0.029      0.02     -0.01    -0.025      0.02     0.005         0         0     -0.01         0    -0.014    -0.019
    0.13     0.034      0.02    -0.014    -0.027     0.024     0.002    -0.003     0.001     -0.01         0    -0.013    -0.022
    0.14      0.04     0.027    -0.013    -0.034     0.026     0.007    -0.006    -0.012    -0.021    -0.004    -0.019    -0.025
    0.15     0.033     0.017    -0.017    -0.025     0.024         0    -0.009    -0.019    -0.033    -0.009    -0.027     -0.02

VN = T1.Properties.VariableNames;
EKG = T1{:,2:end};
Ts = mean(diff(T1.Time))                        % Sampling Interval (s)

Ts = 0.0100

Fs = 1/Ts                                       % Sampling Frequency (Hz)

Fs = 100

Fn = Fs/2                                       % Nyquist Frequency (Hz)

Fn = 50

tilenr = [1:2:12 2:2:12];

figure
tiledlayout(6,2)
for k = 1:size(EKG,2)
    nexttile(tilenr(k))
    plot(T1.Time, EKG(:,k))
    grid
    title(VN{k+1})
end
sgtitle('EKG Original')

[FTs1,Fv] = FFT1(EKG,T1.Time);

figure
tiledlayout(6,2)
for k = 1:size(FTs1,2)
    nexttile(tilenr(k))
    plot(Fv, abs(FTs1(:,k)*2))
    grid
    title(VN{k+1})
end
sgtitle('EKG Original — Fourier Transform')

Considering the leads individually —

for k = 1:size(EKG,2)
    figure
    tiledlayout(1,2)
    nexttile
    plot(T1.Time, EKG(:,k))
    grid
    title(VN{k+1})
    nexttile
    plot(Fv, abs(FTs1(:,k)*2))
    grid
    title("Spectrum "+string(VN{k+1}))
end

% return

EKGfilt = sgolayfilt(EKG, 4, 7);                           % Savitzky-Golay Filter

% figure
% tiledlayout(6,2)
% for k = 1:size(EKG,2)
%     nexttile(tilenr(k))
%     plot(T1.Time, EKGfilt(:,k))
%     grid
%     title(VN{k+1})
% end
% sgtitle('EKG Filtered')
% 
% [FTs2,Fv] = FFT1(EKGfilt,T1.Time);
% 
% figure
% tiledlayout(6,2)
% for k = 1:size(FTs2,2)
%     nexttile(tilenr(k))
%     plot(Fv, abs(FTs2(:,k)*2))
%     grid
%     title(VN{k+1})
% end
% sgtitle('EKG Filtered — Fourier Transform')
% 
% figure
% tiledlayout(2,2)
% nexttile
% plot(T1.Time, T1.II)
% grid
% title('Lead II, Original')
% nexttile
% plot(T1.Time, EKGfilt(:,2))
% grid
% ylim([min(ylim) 1])
% title('Lead II, Filtered')
% nexttile
% plot(Fv, abs(FTs1(:,2))*2)
% grid
% title(["Lead II Original" "Fourier Transform"])
% nexttile
% plot(Fv, abs(FTs2(:,2))*2)
% grid
% title(["Lead II Filtered" "Fourier Transform"])
% sgtitle('Lead II — Detail')


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

The frequency resolution (with a Nyquist frequency of 50 Hz) means that is impossible to specifically detect 50 Hz or 60 Hz mains/power-line noise in these signals. If a proper anti-aliasing filter was not used ahead of the ADC stage, then those frequencies are aliased into the recording and there is no way to remove them at this point. If an appropriate anti-aliasing filter was used (as well as the usual steps taken to prevent those from being recorded, such as a proper reference electrode), then they are simply not in the record and there is no need to remove them.

Most of the leads are relatively ‘clean’, hoiwever there are some exceptions, and I do not see any significant baseline variation in any of them.

Lead III has a significant amount of broadband noise. Wavelet denoising or a Savitzky-Golay filter is the only way to deal with that.

There is some high-frequency noise of an unknown sort in the augmented leads (, and ) however it is relatively low amplitude and I could get the requisite interval information from the other leads, if I was doing a clilnical interpretation of this trace, since the only characteristics I would look for in those leads are the R and T deflections.

Applying a 40 Hz lowpass filtter to them produces —

Leads_4_6_filt = lowpass(EKG(:,1:6), 40, Fs, 'ImpulseResponse','iir');
[FT_Leads_4_6_filt, Fv] = FFT1(Leads_4_6_filt,T1.Time);

for k = 1:size(Leads_4_6_filt,2)
    figure
    tiledlayout(1,2)
    nexttile
    plot(T1.Time, Leads_4_6_filt(:,k))
    grid
    title(VN{k+1})
    nexttile
    plot(Fv, abs(FT_Leads_4_6_filt(:,k)*2))
    grid
    title("Spectrum "+string(VN{k+1}))
end

Applying a Savitzky-Golay filter to them produces —

Leads_4_6_SGfilt = sgolayfilt(EKG(:,1:6), 4, 7);
[FT_Leads_4_6_SGfilt, Fv] = FFT1(Leads_4_6_SGfilt,T1.Time);

for k = 1:size(Leads_4_6_filt,2)
    figure
    tiledlayout(1,2)
    nexttile
    plot(T1.Time, Leads_4_6_SGfilt(:,k))
    grid
    title(VN{k+1})
    nexttile
    plot(Fv, abs(FT_Leads_4_6_SGfilt(:,k)*2))
    grid
    title("Spectrum "+string(VN{k+1}))
end

It is really not possible to do much more than that with them and still retain the essential characteristics necessary to properly interpret that trace.

.

[rawaa mejri]

@Star Strider Thanks a lot,

In the dataset, there are only the four types of noise I mentioned (electrodes_prb, burst_noise, static_noise, and baseline_drift).

During signal extraction, I used this formula: `data(:, i) = (data(:, i) - baseline(i)) / gain(i);` (as we discussed previously).

Regarding high-frequency noise , broadband noise, and others not mentioned in the dataset, have they been addressed? If not, how should I handle them? Specifically for static noise, I'm focusing on the code you provided in your previous comment.

[Star Strider]

As always, my pleasure!

I have done my best to address them. There is no specific way to handle each different type of noise, since they all show up in the signal without any sort of designation. There are limited ways of addressing noise in general, and that is with various types of filtering, or using wavelet denoising.

Using a highpass filter is the best way to eliminate baseline variation. No other baseline processing should be necessary after that, since that will remove constant offset and low-frequency baseline variation.

You can also use wavelet processing to extract various components of the EKG signal. (I have not done that here. I have very little experience with it.)

You can set all of the leads to the same scale using the rescale function, if that is what you want to do. (I would not do that, unless you use the maximum and minimum values in each lead, because the relative amplitudes of the various deflections in the different leads is clinicallly important.)

[rawaa mejri]

Thank you very much @Star Strider.

Do you advise me to treat each type of noise separately or to use the same treatment you did here and ignore the noise in the dataset? After that, should I apply Z-score normalization?

Thank you very much.

[Star Strider]

As always, my pleasure!

I doubt that there is a way to treat each sort of noise differently, essentially because there is no way to differentiate them. You could see if using wavelets to detect them individually would work, however thaty is likely not possible either.

I am not certain why you would want to normalise them at all. The impiortant features in an EKG (and for most — if not all — physiological signals) is the signal morphology (important signal features) in the signals, and normalising them eliminates those features. Also, the relative amplitudes of the R-deflection (and other deflections) in the various leads are important, for example in calculating the cardiac axis, and normalising the various leads to ome set of amplitudes destroys that. There is also normal amplitude variations in any given lead that is also important and shoold not be changed by normalising them.

I would just elimiinate as much high-frequency noise and low-frequency baseline variation and offset as possible, and leave it at that. There are various ways to do that, and I outlined frequency-selective filters and multiband filters (Savitzky-Golay filters) in my discussion here, that are the only available ways to deal with the various types of noise. If you have the Wavelet Toolbox, you can see if wavelet denoising (the wdenoise function, and related functions) can improve these signals, however I doubt if they could significantly improve what has already been done here.

[rawaa mejri]

Hello @Star Strider , I hope your are doing well ,

I use wavelet transform , I test differents wavelets as shown with picture below.

I see that bior3.9 with adaptative thershold is the best ? Do you confirm that ?

Thanks a lot :)

[Star Strider]

I do not have an extensive background with wavelets, although I have used them. I am not familiar with bior splines. If it gives you acceptable results, and has a theoretical basis that makes it appropriate for EKG analysis, then use it! Congratulations on discovering thtat approach.