What filtration should be used for a respiratory signal between 5 and 60 breaths per minute?

I have applied LPF and HPF filtering, but for a low-frequency respiratory signal, the signal after HPF is distorted and edge detection is incorrect, as can be seen in the attached photo. Is it possible to use some kind of filtration that can handle the respiratory signal range of 5-60 breaths per minute? Is it better to use a findpeaks function with appropriate limitations and only LPF filtering as in diagram 2 in the attached photo?

% ------------- LPF ------------------------------
N     = 5;     % Order
Fstop = 1.4;  % Stopband Frequency 
Astop = 30;    % Stopband Attenuation (dB)
Fs    = 25;  % Sampling Frequency
h = fdesign.lowpass('n,fst,ast', N, Fstop, Astop, Fs);
hfiltLP = design(h, 'cheby2', 'SystemObject', true);


% ------------- HPF ------------------------------
N     = 4;     % Order
Fstop = 4/60;     % Stopband Frequency min 4 oddechow na minute
Astop = 60;    % Stopband Attenuation (dB)
Fs    = 25;  % Sampling Frequency
h = fdesign.highpass('n,fst,ast', N, Fstop, Astop, Fs);
hfiltHP = design(h, 'cheby2', 'SystemObject', true);

[Nader]

% -------------- LPF ------------------------------ N = 5; Fstop = 1.4; Astop = 30; Fs = 25; h = fdesign.lowpass('n,fst,ast', N, Fstop, Astop, Fs); hfiltLP = design(h, 'cheby2', 'SystemObject', true);

% -------------- HPF ------------------------------ N = 4; Fstop = 4/60; Astop = 60; Fs = 25; h = fdesign.highpass('n,fst,ast', N, Fstop, Astop, Fs); hfiltHP = design(h, 'cheby2', 'SystemObject', true);

Answer 1

[William Rose]

@Pawel Slazak,

%% Load data
a=importdata('RAW_signal.txt');
x=a.data(:,1);              %x=channel 1
fs=25;                      %sampling frequency (Hz)
t=(0:length(x)-1)/fs;       %time vector
%% Filter data
fco=1.4;                    %lowpass cutoff (Hz)
[b,a]=butter(4,fco/(fs/2)); %4th order Butterworth lowpass
y=filter(b,a,x-mean(x));    %filter the data
y=y+mean(x);                %add back the mean value
%% Find minima
[pks,locs]=findpeaks(-y,fs,'MinPeakProminence',10);  %find peaks and peak locations
%% Compute rate
rRate=60./diff(locs);       %rate associated with each breath (1/minute)
tRate=locs(2:end);          %ending time of each breath
%% Plot results
figure; subplot(211)
plot(t,x,'-r',t,y,'-b',locs,-pks,'k^'); 
title('Respiration'); grid on; xlabel('Time (s)'); ylabel('Amplitude (mV)');
legend('Raw','Filtered','Minima','Location','Southeast');
subplot(212); stairs(tRate,rRate,'-r'); grid on; 
title('Respiratory Rate'); xlabel('Time (s)'); ylabel('Rate (bpm)');

The code above finds and plots the the minima. It also computes and plots the respiratory rate versus time. I chose to associate each breath's rate with the ending time of each breath, and I used the stairs() plot for the rate, because that is how a real time app would work. You said you hope to implement this in real time.

[Pawel Slazak]

Thanks for your help.

Answer 2

[William Rose]

@Pawel Slazak, pease attach the raw data.

Please explain what your goal is, because this will determine what approach to take.

If it is your goal to determine the respiratory rate, then it does not matter if we use upward or downward slope of the signal to mark each breath, as long as we find a consistent way that is not affected by the noise on the signal.

If it is your goal to identify the time of the start of each inspiration, then please specify whether more positive values of the signal correspond to greater or lesser lung volume. We will want to identify the minimum points of the lung volume signal.

It will not be hard to do this. We just need to choose the right filter setting, and we must decide if we will use peaks in filtered dV/dt to do it - which depends on your answers to the quesitons above.

Please upload the raw data, as a text file or as a .mat file.

[Pawel Slazak]

My goal is to estimate the respiratory rate in a real-time system.

In the attached .txt file there are two signals, "Sensor_1" and "Sensor_2". Sensor 2 was incorrectly mounted, so it is not taken into consideration.

[William Rose]

@Pawel Slazak, your plots suggest that you may be able to estimate the respiratory rate effectively by finding the minima of the low-pass filtered signal. Do this by finding the peaks of the negative of the low pass filtered signal. You may also want to use the MinPeakHeight criterion or the MinPeakProminence criterion, to avoid unwanted minima.

[William Rose]

@Pawel Slazak,

Thank you. Can you upload a file of raw data, such as the raw data for the upper plot in the figure you posted?

[Pawel Slazak]

This is the signal I posted, but in the attached figure I have included only part of the signal due to its size.

[William Rose]

Thank you for posting the data. I will look at it later tonight.

[William Rose]

@Pawel Slazak,

I will ignore the channel 2 data since it appears to be quite a sdifferent signal.

Channel 1 is the raw data from your original plot, but for a longer diraiton. The period is about 12.5 seconds initially and 1.0 seconds at the end of the 640 second recording. This corresponds 5 to 60 breaths per minute. This is a very wide a range of respiratory rates - beyond the normal physiological range. What animal is this? What are the experimental conditions?

[Pawel Slazak]

The signal was obtained artificially in order to write an algorithm.

[William Rose]

@Pawel Slazak, Thanks. That makes sense.