How to add variable High Frequency Noise and Baseline drift to signals in a dataset

I have a dataset of ecg signals. I have load it all of them into the workspace and filtered with a band-pass butterworth filter. Now I have to add high frequency noise and baseline drift to these filtered signals. The total noise added (high frequency noise + drift) needs to be variable and not be the same.

For the noise can I use a Gaussian or in this case it is wrong? I tried to use awgn, but if it wrong in this context, tell me.

How can I do it?

This is my code:

%% load data and filtering the signals
clear; clc; % clear workspace
Fs = 360; % [Hz]
time_interval = 10; % [sec]

files = dir(fullfile("dataset/","*.mat")); % all dataset files
[b,a] = butter(3,[1 30]/(Fs/2),"bandpass"); % bandpass butterworth filter 
filt_sn = zeros(numel(files), 3600); % prealloc 

%--load all the data and filter it--%
for i = 1:numel(files)
    data(i) = load(fullfile("dataset/",files(i).name)); % load all data
    filt_sn(i,:) = filtfilt(b,a,data(i).val); % filter all signals
end

%--Adding Noise and Baseline Drift--%
% prealloc for noise %
noisy_sn = zeros(numel(files), 3600); % prealloc 
SNR = zeros(1,numel(files)); %prealloc
% prealloc for drift
dritf_sn = zeros(numel(files), 3600); % prealloc
A = zeros(1,numel(files)); % prealloc
N = zeros(1,numel(files)); % prealloc
minDriftOffset = zeros(1,numel(files)); % prealloc
maxDriftOffset = zeros(1,numel(files)); % prealloc

%--MY ATTEMPT SO FAR--%
for j = 1:numel(files)
    %--add noise--%
    SNR(j) = 8 + (15 - 8).*rand(1);
    noisy_sn(j,:) = awgn(filt_sn(j,:),SNR(j),'measured');
    %--add baseline drift--%
    %--generate random variable for all signals--%
    X = linspace(0,2*pi,length(filt_sn(j,:))); 
    A(j) = -200 + (200 - (-200)).*rand(1);
    N(j) = 50 + (10 - 50).*rand(1);
    minDriftOffset(j) = -100 + (0 - (-100)).*rand(1);
    maxDriftOffset(j) = 10 + (150 - 10).*rand(1);
    %--sum up to distort the signal--%
    cos_wave = A(j)*cos(X(j)./N(j));
    slantedLine = linspace(minDriftOffset(j),maxDriftOffset(j),length(filt_sn(j,:)));
    dritf_sn(j,:) = noisy_sn(j,:)+ cos_wave + slantedLine;
end

%--plot the signal/s you want--%
%--to plot: chanhe index 'i' in (i,:);
% figure(); subplot(2,1,1); plot(T,noisy_sn(1,:)); grid on;
% title("Noisy ECG signal"); xlabel("time [sec]"); ylabel("voltage [mV]");
%
% figure(); subplot(2,1,1); plot(T,drfit_sn(1,:)); grid on;
% title("Distorted ECG signal (noise + drift)"); xlabel("time [sec]"); ylabel("voltage [mV]");

[Eamon Gekakis]

could injecting noise with a sin curve + specified high frequency sample rate + randi() generated noise suffice? Drift can be added by multipling your data by a deviation slope curve.

Answer 1

[Image Analyst]

There are lots of kinds of noise. Theoretically what kind does your signal typically have? You can use rand() and randn(). You can also change the amplitude and range of the noise on an element-by-element basis if you want. Do you want to do that, or is your noise relatively stationary across your signal?

For the baseline moving, you have to decide how you want the baseline to move. Do you just want a linear ramp added to the signal? Do you want a sinusoid added to the signal? Do you want the baseline to move fairly randomly, like 1-D Perlin noise was applied to it?

Google Search Perlin Noise 1D

[Sabaudian]

I have to add a high frequency noise (this is my assignment, I don't have a specific frequency or other data) to all the ECG signals and this noise needs to be different for every signals (noise1 for signal1, noise2 for signal2,... I have 48 signals in my dataset).

For the noise, I tried with awgn function, but I don't know if "my noise" is a high frequency noise.

For the baseline drift, I think I have solved it, as you can see in my attempt above.

[Image Analyst]

I think using rand will give you all frequencies of noise. To get only the high frequencies you'd have to do an FFT and zero out the elements below some specified element/frequency.

[Sabaudian]

ok, but how can I "zero out the elements below some frequency" like you said?

[Sabaudian]

I made another attempt on a single ecg signal to testing, but I really don't understand your prievious suggestion

My attempt:

clear; clc; % clearing the workspace
%--Load Real signal--%
load ('100m.mat');
ecg = val/200;
Fs = 360; % sample frequency [Hz]
L = length(ecg); % Signal Length
T = linspace(0,L/Fs,L); % time axis
figure (1); plot(T, ecg); grid on; % plot og. signal
title('ECG Signal'); xlabel('time (sec)'); ylabel('voltage (mV)');

%--ATTEMPT--%
fn = 60; % Frequency to add
yn = 0.2*sin(2*pi*fn*T); % noise
y = fft(ecg) + yn; % add noise
noisy_sn = abs(ifft(y)); % final result
figure (2); plot(T, noisy_sn); grid on; % plot final result
title('Noisy ECG Signal'); xlabel('time (sec)'); ylabel('voltage (mV)');

[Image Analyst]

That's basically adding 60 Hz "power line" noise to the fft spectrum (instead of the time domain signal) which doesn't make sense. Why would a pure deterministic signal like that exist in the frequency domain? If you just wanted to add line noise to the signal you'd either add that since wave to the time domain signal, or you'd make a "spike" at the one single element that represented 60 Hz.

But maybe if you want high frequency noise more than 100 Hz, you can try this:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;
markerSize = 30;
fprintf('Beginning to run %s.m ...\n', mfilename);
%--Load Real signal--%
load ('100m.mat');
ecg = val/200;
Fs = 360; % sample frequency [Hz]
L = length(ecg); % Signal Length
T = linspace(0,L/Fs,L); % time axis
F = linspace(-Fs/2, Fs/2, L); % Frequency axis
subplot(3, 2, 1);
plot(T, ecg, 'b-'); 
grid on;
title('ECG Signal', 'FontSize', fontSize); 
xlabel('time (sec)', 'FontSize', fontSize); 
ylabel('voltage (mV)', 'FontSize', fontSize);
% Save y axis limits for later.
yl = ylim;
% Get the spectrum of the original signal.
ftEcg = fft(ecg);
shiftedFT = fftshift(abs(ftEcg));
subplot(3, 2, 2);
plot(F, shiftedFT, 'b-'); 
grid on;
title('Spectrum of original ECG Signal', 'FontSize', fontSize); 
xlabel('Frequency (Hz)', 'FontSize', fontSize); 
ylabel('Signal', 'FontSize', fontSize);
% Make noise vector
amplitude = 2;
noiseSpectrum = amplitude * rand(1, length(T));
shiftedNS = fftshift(noiseSpectrum);
% Zero out less than 100 Hz.
index = round(length(shiftedNS) * 100 / Fs)
middleIndex = length(shiftedNS) / 2
shiftedNS((middleIndex - index) : (middleIndex + index)) = 0;
subplot(3, 2, 3);
plot(F, shiftedNS, 'b-'); 
grid on;
title('Spectrum of Noise Signal', 'FontSize', fontSize); 
xlabel('Frequency (Hz)', 'FontSize', fontSize); 
ylabel('Signal', 'FontSize', fontSize);
% Add noise spectrum to original signal spectrum
noisySpectrum = abs(shiftedFT + shiftedNS); 
subplot(3, 2, 4);
plot(T, noisySpectrum, 'b-');
grid on;
title('Noisy ECG Spectrum', 'FontSize', fontSize); 
xlabel('Frequency (Hz)', 'FontSize', fontSize); 
ylabel('signal', 'FontSize', fontSize);
% Shift it back to get the spectrum in a form where we can use ifft:
noisySpectrum = ifftshift(noisySpectrum);
subplot(3, 2, 5);
plot(T, noisySpectrum, 'b-');
grid on;
title('Noisy ECG Spectrum', 'FontSize', fontSize); 
xlabel('Frequency (Hz)', 'FontSize', fontSize); 
ylabel('signal', 'FontSize', fontSize);
% Return to time domain
noisyTimeDomainSignal = ifft(noisySpectrum);
subplot(3, 2, 5:6);
plot(T, noisyTimeDomainSignal, 'b-');
grid on;
title('Noisy ECG Signal', 'FontSize', fontSize); 
xlabel('time (sec)', 'FontSize', fontSize); 
ylabel('voltage (mV)', 'FontSize', fontSize);
ylim(yl)

No guarantees - you'll have to review my math.