ECG time domain synchronized averaging

[bababui]

Hello,

I am trying to understand the synchronized averaging technique and implement it in MATLAB. I have a dataset which represents a noisy ECG signal with 10 epochs of the signal, and trying to apply the above mentioned method to filter it.

In general I know that I'm supposed to separate the given ECG signal into epochs line them up and add them together and then divide by the number of epochs to get a cleaner epoch.

Can someone provide some detail, I think I'm missing some important points.

 

[werewolf]

I think what you have described is the traditional averaging process used in signal processing and image processing to reduce noise which is different than (TSA) time domain synchronized averaging.

In the traditional averaging you have a time domain signal with noise which was measured using some digital analyzer.
e.g.: [ 2 -3 5 -1 4 6 -7 5 -9 10 0 -3……………….]
Where every number in this signal represents a measurement made in one sample (I used exact integers so it is easier to do calculations but in general the signal will be real numbers).

In the traditional averaging:
1-Choose n (odd number that represent the size of the averaged block)
2-Padding the signal (add zeros to the start and the end of the signal e.g.: 0 0 0 2 -3 5 -1 ...0 0 0)
3-Replace every measurement by its average with its neighborhood.(you can always remove the padding after you finish your calculations!)
e.g.:
Assume n=3 (so we will look for one reading before and one reading after[moving average])
Data= [ 0 2 -3 5 -1 4 6 -7 5 -9 10 0 -3……………….]
A_Data= [ 0 -1 4 1 8 9 3 4 ………………….……………….]/3
But for the TSA we have to have two signals, the first one is your data readings and the other one is the trigger.
-The reason for this is based on the idea of cause and effect.
e.g.: a machine gives you a loud noise so it has a problem and you couldn’t identify where is the problem using the noise readings only so you measure two things at once the noise readings and the rpm of the part which you doubt to be the source of the noise.
Then you apply the TSA to reduce the noise if and only if this noise pattern is related to the rpm of this part…

For more information about TSA I recommend searching in Google

 

[ramyapc]

Hi, can you please tell me hw can i average ecg signal. i have a signal with ten beats, all these beat is to be averaged and a signal beat is needed.

 

[werewolf]

time domain synchronized averaging is the most suitable for this kind of signal(I mean of course the ECG signal)

but to do so we have to create our trigger signal which we can assume simply to be zeros every where and one at the place of the peak(the top part of the R signal)

Now because we have the two signals we can use TSA Easily by simply divide the beats signal into periods of time corresponding to the trigger signal.

the next step is normalizing the small signals resulting from the first step so that they all have the same time period(compress long signal and stretch short signal so that all the signals has the same time frame)

now for the averaging step simply add the signals by adding every point to the points corresponding to the same time in the other beats signals and finally dividing by the number of the small beats signals you have your time domain averaged signal...

of course you can try low-pass filters on this signal but it will need more experience in filtering than using TSA.

Simple algorithm:
1-your signal will be y(t) where t is the time vector and y is the ECG reading
2-your trigger will be pulse(t) which is zeros every where except for the maximum R peak.
3-the number of beats at your signal is n
4-use the trigger to divide the time vector into t1, t2,.......tn vectors
5-use the trigger to divide y(t) to y1(t1), y2(t2),.........yn(tn) signals.
6-Normalize the signal in the time axis direction so that all the signals will be represented in the same period tp: y1(tp), y2(tp),.....yn(tp)
7-your averaged signal will equal y_av(tp)=(y1(tp)+y2(tp)+......+yn(tp))/n
8-Now you have your averaged signal that is more smother than any of the beats and also more representative.

---------- Post added at 21:20 ---------- Previous post was at 21:16 ----------

Those are also images that may help visualize the algorithm while reading the steps:

**broken link removed**

Also I must say that there is a lot more in ECG filtering than just that but I'm trying to simplify things.
I'm using TSA now in vibration analysis which clearly states that it's a general method and can be used in many fields and off course there might be a methods that is specialized in ECG filtering that I don't know about or just know the main Idea with no experince

 

[werewolf]

This MatLab function may work just fine to do time domain synchronized averaging,

for this function to work you need to prepare three vectors
-time vector representing the time samples.
-pulse vector representing the dividing of your signal into beats(zeros every where and 1 at the start of the beat)
-and finally your signal

(they all must have the same dimensions)

%Code



function [] = TSA(time,pulse,y)
%This function calculates the time domain averaged of a signal
% time = your time vector;
% pulse= your trigger signal;
% y = your signal;

% The number of pulses in one cycle
ppc=1;
% Locations of positive spikes, vector containing the nr:s of the samples where spikes are located
pulse_positions = find(pulse == 1);
% Number of samples in each gap between two pulses
samp_per_gap = diff(pulse_positions);
% Maximum nr of samples in the longest gap
max_samp_per_gap = max(samp_per_gap);
% Number of gaps
nr_gap = length(samp_per_gap);
% Number of complete revolutions, 2 pulses per revolution
nr_rev = fix(nr_gap/ppc);
% Maximum nr of samples in one revolution (only upsampling allowed)
samp_per_rev = ppc*max_samp_per_gap;
% Distance between equidistant angle samples, 360 degreed per revolution
delta_angle = 360/samp_per_rev;
% New equidistant angle vector
angle = (0:nr_rev*samp_per_rev-1)*delta_angle;
% Vector with angle location of spikes
angle_spikes = angle(1:samp_per_rev/ppc:end);
% Truncate pulse_positions to fit angles_spikes, keep only complete revolutions
pulse_positions = pulse_positions(1:length(angle_spikes));
% Interpolate to find anglular position at each time sample
angles_of_t_sampl = interp1(time(pulse_positions), angle_spikes, time, 'cubic', 'extrap');
% Interpolate to find signal values at equidistant angle samples
acc0_resampled = interp1(angles_of_t_sampl, pulse, angle);
acc1_resampled = interp1(angles_of_t_sampl, y, angle);

% Calculate synchronised time domain average

% Angle vector for one revoluion
angle_rev = angle(1:samp_per_rev);
% Define acc_mean vector
acc0_mean = zeros(1,samp_per_rev);
acc1_mean = zeros(1,samp_per_rev);
% Loop to calculate the synchronised time domain average
for M = 1:nr_rev; acc0_mean = acc0_mean + acc0_resampled(1+(M-1)*samp_per_rev:M*samp_per_rev); end
acc0_mean = acc0_mean/nr_rev;
for M = 1:nr_rev; acc1_mean = acc1_mean + acc1_resampled(1+(M-1)*samp_per_rev:M*samp_per_rev); end
acc1_mean = acc1_mean/nr_rev;

% Plot synchronised time domain average
figure(1);
subplot(3,1,1);plot(angle_rev,acc0_mean),title('Pulse averaged time domain signal');
subplot(3,1,2);plot(angle_rev,acc1_mean),title('Averaged time domain signal');

 

[albara]

This MatLab function may work just fine to do time domain synchronized averaging,

for this function to work you need to prepare three vectors
-time vector representing the time samples.
-pulse vector representing the dividing of your signal into beats(zeros every where and 1 at the start of the beat)
-and finally your signal

(they all must have the same dimensions)

%Code



function [] = TSA(time,pulse,y)
%This function calculates the time domain averaged of a signal
% time = your time vector;
% pulse= your trigger signal;
% y = your signal;

% The number of pulses in one cycle
ppc=1;
% Locations of positive spikes, vector containing the nr:s of the samples where spikes are located
pulse_positions = find(pulse == 1);
% Number of samples in each gap between two pulses
samp_per_gap = diff(pulse_positions);
% Maximum nr of samples in the longest gap
max_samp_per_gap = max(samp_per_gap);
% Number of gaps
nr_gap = length(samp_per_gap);
% Number of complete revolutions, 2 pulses per revolution
nr_rev = fix(nr_gap/ppc);
% Maximum nr of samples in one revolution (only upsampling allowed)
samp_per_rev = ppc*max_samp_per_gap;
% Distance between equidistant angle samples, 360 degreed per revolution
delta_angle = 360/samp_per_rev;
% New equidistant angle vector
angle = (0:nr_rev*samp_per_rev-1)*delta_angle;
% Vector with angle location of spikes
angle_spikes = angle(1:samp_per_rev/ppc:end);
% Truncate pulse_positions to fit angles_spikes, keep only complete revolutions
pulse_positions = pulse_positions(1:length(angle_spikes));
% Interpolate to find anglular position at each time sample
angles_of_t_sampl = interp1(time(pulse_positions), angle_spikes, time, 'cubic', 'extrap');
% Interpolate to find signal values at equidistant angle samples
acc0_resampled = interp1(angles_of_t_sampl, pulse, angle);
acc1_resampled = interp1(angles_of_t_sampl, y, angle);

% Calculate synchronised time domain average

% Angle vector for one revoluion
angle_rev = angle(1:samp_per_rev);
% Define acc_mean vector
acc0_mean = zeros(1,samp_per_rev);
acc1_mean = zeros(1,samp_per_rev);
% Loop to calculate the synchronised time domain average
for M = 1:nr_rev; acc0_mean = acc0_mean + acc0_resampled(1+(M-1)*samp_per_rev:M*samp_per_rev); end
acc0_mean = acc0_mean/nr_rev;
for M = 1:nr_rev; acc1_mean = acc1_mean + acc1_resampled(1+(M-1)*samp_per_rev:M*samp_per_rev); end
acc1_mean = acc1_mean/nr_rev;

% Plot synchronised time domain average
figure(1);
subplot(3,1,1);plot(angle_rev,acc0_mean),title('Pulse averaged time domain signal');
subplot(3,1,2);plot(angle_rev,acc1_mean),title('Averaged time domain signal');

Thanks for your code but I could not use it for my signal and I do not know how to apply it on different signal. In addition to that, I do not know what you mean by trigger signal

 

[albara]

Thanks for your code but I could not use it for my signal and I do not know how to apply it on different signal. In addition to that, I do not know what you mean by trigger signal