Extracting individual breathing cycles from a signal and plotting the mean cycle

I have a breathing trace, which I have smoothed and use FindPeaks to identify the peaks. How do I extract the individual beathing cyles between the peaks so that I can plot these over each other and calculate the mean breathing cycle?

I know there is a baseline shift so I also need to look at detrending the data first.

I would be very grateful for any input or direction to where I might find information on how to do this.

I have tried to extract the first 3 regions of interest but an 'unrecognized function error'.

roilims = [0 227; 227 480; 480 745];

>> sigroi = extractsigroi(dt_B,roilims);

Unrecognized function or variable 'extractsigroi'.

Answer 1

[Star Strider]

If the peaks define the cycle, then finding the ‘valleys’ between them will define the beginning and end points of each cycle. I would use either findpeaks with the negative of the signal (turning the ‘valleys’ into peaks) or the islocalmin function.

An example of that to create an ensemble average (not necessary, use only the parts of the code that apply to your problem) is Combining repetitive curves into one average curve and specifically this Comment.

.

[Gail Distefano]

Thanks very much. This is really helpful. I have given it a try. I attach the a file with my trace, sampled at a frequency of 30Hz.

I have run your code and get these images. Please may you explain the significance of '-55' and '+320'. I seem to be extracting more than one cycle.

'if (plocs(k)+320) <= numel(t)

idxrng = max(1,plocs(k)-55):min(numel(t),plocs(k)+320);'

[Star Strider]

My pleasure!

I can help to implement it with your data if you care to share your data.

EDIT — (24 Feb 2022 at 19:52)

Your data are different from the previous data, so it was necessary to tweak the code a bit.

LD = load('breath_smoothed.mat');
dt_B = LD.dt_B;
y = dt_B;
Fs = 30;
t = linspace(0, numel(y)-1, numel(y))/Fs;                           % Assume 1kHz Sampling Frequency (Supply Correct Time Vector)
Ts = mean(diff(t));

FTy = fft(y-mean(y))/numel(y);                                      % Calculate & Plot Fourier Transform
Fv = linspace(0, 1, numel(t)/2+1)*Fs/2;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTy(Iv))*2)
grid
xlim([0 1.5])
xlabel('Frequency')
ylabel('Amplitude')

yf = highpass(y, 3E-3*Fs, Fs, 'ImpulseResponse','iir');                                     % Eliminate Baseline Variation

[pks,plocs] = findpeaks(yf, 'MinPeakProminence',0.5);               % Peaks & Location Indices
[vls,vlocs] = findpeaks(-yf, 'MinPeakProminence',0.5);              % Valleys & Location Indices

figure
plot(t, yf)                                                         % Detail Of Original Waveform
hold on
plot(t(plocs),pks,'^r') 
hold off
% xlim([0 2.5])
xlabel('Time (s)')
ylabel('Amplitude (Unit)')

Vidx = mean(diff(vlocs));
for k = 1:numel(vls)-1
    idxrng{k,:} = max(1,vlocs(k)):min(numel(t),vlocs(k+1));         % Generate Ensemble From Individual Records
    ensb{k,:}(1:numel(idxrng{k})) = yf(idxrng{k});
end

idxrng_max = max(cellfun(@numel,idxrng));                           % Generate MAtrices For Statistics Plots
idxrng_mtx = zeros(size(idxrng,1), idxrng_max);
ensb_mtx = idxrng_mtx;

figure
hold on
for k = 1:size(ensb,1)                                              % Plot Data, Populate Matrices
    plot(idxrng{k}, ensb{k})
    idxrng_mtx(k,1:numel(idxrng{k})) = idxrng{k};
    ensb_mtx(k,1:numel(ensb{k})) = ensb{k};
end
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude (Unit)')

figure
hold on
for k = 1:size(ensb_mtx,1)                                          % Plot Ensemble
    plot(1:numel(idxrng{k}), ensb{k})
end
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude (Unit)')

ensb_mean = mean(ensb_mtx);
ensb_std = std(ensb_mtx);
tci = tinv(0.975,size(ensb,1)-1);                                   % 95% Critical t-Statistic
err = tci*ensb_std/sqrt(size(ensb,1));                              % 95% Confidence Intervals

figure
plot(t(1:numel(ensb_mean)), ensb_mean, '-b')
hold on
plot(t(1:numel(ensb_mean)), ensb_mean+err, '--r')
plot(t(1:numel(ensb_mean)), ensb_mean-err, '--r')
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude (Unit)')
legend('Ensemble Average', '±95% Confidence Intervals')

Figure 5 —

It can be tweaked further if necessary if you have a different desired result.

.

[Gail Distefano]

Hi,

I thought I attached my data file above. Is it something else you need? It is also important for me to understand the script too because I going to be acquiring lots of breathing traces for different patients and will need to be able to adapt the script.

'breath_smoothed.mat'

Thanks very much again.

[Gail Distefano]

Thanks very much again for your prompt reply!

[Star Strider]

My pleasure!

It took a bit of tweaking to get my previous code to work with your data.

If my Answer helped you solve your problem, please Accept it!

.

[Gail Distefano]

Hi, I have been clinical the last few days so only just finding time to try the code. I am getting an error half way through. I runs well till this point, would be very gratefu if you can still help. I thought maybe there is an error in the second line:for k = 1:numel(vls)-1. 'vls' should be 'vlocs'? But I still get the same error.

Vidx = mean(diff(vlocs));

for k = 1:numel(vls)-1

idxrng{k,:} = max(1,vlocs(k)):min(numel(t),vlocs(k+1)); % Generate Ensemble From Individual Records

ensb{k,:}(1:numel(idxrng{k})) = yf(idxrng{k});

end

Unable to perform assignment because brace indexing is not supported for variables of this type.

[Star Strider]

In my code, ‘vls’ (valley magnitudes) and ‘vlocs’ (the index locations of them) should be the same size, so that is not the problem.

With respect to this error:

Unable to perform assignment because brace indexing is not supported for variables of this type.

My code runs without error with the .mat file provided. I have no idea what is throwing that error.

If you are using a different .mat file, post it and I will see if I can figure out what the problem is with it.

[Gail Distefano]

Many thanks again, I cleared the workspace and started from scratch and it works well. Apologies!

I have tried to run it with another dataset at 25hz. I attach the file here. It runs fine until Figure 3, when it plots only part of the trace. I did change fs=25. It would be really helpful if you are able to tell me why so that I can learn to tweak the code myslef. Thanks very much

[Star Strider]

It took me a while to figure out what the problem was that you referred to.

My code plots the identified waves defined as those that findpeaks returns both ‘valley’ indices for in ‘vlocs’. If a peak (or valley) cannot be identified as such since it is at the end of the vector and does not have a defined peak (or valley), findpeaks will not detect it. In this instance, the last peak does not have an identified trailing ‘valley’. One way to deal with that is to append the last index of the vector to the ‘vlocs’ vector:

[vls,vlocs] = findpeaks(-yf, 'MinPeakProminence',0.5); % Valleys & Location Indices

vlocs = [vlocs, numel(yf)];

The loop that then detects and returns the respiratory cycles is then modified to be:

Vidx = mean(diff(vlocs));

for k = 1:numel(vlocs)-1

idxrng{k,:} = max(1,vlocs(k)):min(numel(t),vlocs(k+1)); % Generate Ensemble From Individual Records

ensb{k,:}(1:numel(idxrng{k})) = yf(idxrng{k});

end

It does not make a significant difference in the result, however it now includes even the incomplete respiratory cycles.

The full code for this file is now:

LD = load('BreathData_Type2.mat');

y = LD.C;

Fs = 25;

t = linspace(0, numel(y)-1, numel(y))/Fs; % Assume 1kHz Sampling Frequency (Supply Correct Time Vector)

Ts = mean(diff(t));

FTy = fft(y-mean(y))/numel(y); % Calculate & Plot Fourier Transform

Fv = linspace(0, 1, numel(t)/2+1)*Fs/2;

Iv = 1:numel(Fv);

figure(1)

plot(Fv, abs(FTy(Iv))*2)

grid

xlim([0 1.5])

xlabel('Frequency')

ylabel('Amplitude')

yf = highpass(y, 3E-3*Fs, Fs, 'ImpulseResponse','iir'); % Eliminate Baseline Variation

[pks,plocs] = findpeaks(yf, 'MinPeakProminence',0.5); % Peaks & Location Indices

[vls,vlocs] = findpeaks(-yf, 'MinPeakProminence',0.5); % Valleys & Location Indices

vlocs = [vlocs, numel(yf)];

figure(2)

plot(t, yf) % Detail Of Original Waveform

hold on

plot(t(plocs),pks,'^r')

hold off

% xlim([0 2.5])

xlabel('Time (s)')

ylabel('Amplitude (Unit)')

Vidx = mean(diff(vlocs));

for k = 1:numel(vlocs)-1

idxrng{k,:} = max(1,vlocs(k)):min(numel(t),vlocs(k+1)); % Generate Ensemble From Individual Records

ensb{k,:}(1:numel(idxrng{k})) = yf(idxrng{k});

end

idxrng_max = max(cellfun(@numel,idxrng)); % Generate MAtrices For Statistics Plots

idxrng_mtx = zeros(size(idxrng,1), idxrng_max);

ensb_mtx = idxrng_mtx;

figure(3)

hold on

for k = 1:size(ensb,1) % Plot Data, Populate Matrices

plot(t(idxrng{k}), ensb{k})

idxrng_mtx(k,1:numel(idxrng{k})) = idxrng{k};

ensb_mtx(k,1:numel(ensb{k})) = ensb{k};

end

hold off

grid

xlabel('Time (s)')

ylabel('Amplitude (Unit)')

figure(4)

hold on

for k = 1:size(ensb_mtx,1) % Plot Ensemble

plot(1:numel(idxrng{k}), ensb{k})

end

hold off

grid

xlabel('Time (s)')

ylabel('Amplitude (Unit)')

ensb_mean = mean(ensb_mtx);

ensb_std = std(ensb_mtx);

tci = tinv(0.975,size(ensb,1)-1); % 95% Critical t-Statistic

err = tci*ensb_std/sqrt(size(ensb,1)); % 95% Confidence Intervals

figure(5)

plot(t(1:numel(ensb_mean)), ensb_mean, '-b')

hold on

plot(t(1:numel(ensb_mean)), ensb_mean+err, '--r')

plot(t(1:numel(ensb_mean)), ensb_mean-err, '--r')

hold off

grid

xlabel('Time (s)')

ylabel('Amplitude (Unit)')

legend('Ensemble Average', '±95% Confidence Intervals')

I am not certain what to suggest for a general solution, other than to not include the incomplete respiratory cycles at eiher end, as the original code does not. Including them as this version does, risks corrupting the ensemble data with incomplete cycles, so I advise against it and advise instead using the original code and not this modified version.

I did not also include the data that began with the first element of the vector for this reason, although to be consistent, it would be necessary to include the first and last indices of the vector in ‘vlocs’. That would create ‘vlocs’ as:

vlocs = [1, vlocs, numel(yf)];

However they are your data and your results, so it is your choice as to how to work with the incomplete cycles.

.

[Gail Distefano]

Hi, I just wanted to say thanks very much for your reply and times taken to assist me. For a strange reason I did not get a notification so have only just seen this reply. It is extremely helpful. Thanks once again

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

Answer 2

[Cris LaPierre]

extractsigroi is part of the Signal Processing Toolbox. It would appear you do not have it installed. If you have access to it through your license, you can add it using the Add-Ons explorer.

[Gail Distefano]

Thank you. But I definitely have Signal Processing Toolbox installed. Could it be that it is not included in Vn2019b?

[Cris LaPierre]

Yes, that is the case. This function was introduced in R2020b.