# How to calculate respiratory rate from a Doppler radar signal

21 views (last 30 days)
Mohan kand on 21 Mar 2023
Edited: Mohan kand on 19 Apr 2023
I have a respiration signal from Doppler radar (see the radar_signal.mat and ). The sampling frequency is 2 KHz, Pulse repetition time is 0.0005 sec. I have no idea what kind of filter I need to apply to detect the respiratory signal. My goal is to use those signals to detect the Respiratory Rate, RR. Below is my MATLAB code I used to calculate the FFT and RR. I modified the code based on the work published in : https://www.sciencedirect.com/science/article/pii/S2352340921009999. Please let me know if I am correct or not.
[file,path] = uigetfile('*.txt');
Data = textread(file, '%s', 'delimiter', '\n');
% Make Data cells empty if numerical
numericalArray = cellfun(@(s) sscanf(s,'%f').' ,Data, 'un', 0);
Fs = 2000; % Sampling Frequency (Hz)
Fn = Fs/2; % Nyquist Frequency
Ts = 1/Fs; % Sampling Time (sec), pulse repitation time
t = linspace(0, L, L)*Ts; % Time Vector (sec)
figure(1)
grid
xlabel('Time')
ylabel('Amplitude')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Signal centralization
% filter
Passband_Frequency_kHz = 2;
% FFT
[radarSpectrumForResp_filtered, ~] = FFT(Filtered, Fs) ;
% Estimate respiratory rate
figure(1)
subplot(3,3,[1 2])
xlabel('time (s)')
ylabel('amplitude (V)')
subplot(3,3,[4 5])
plot(t, Filtered)
ylabel('amplitude (V)')
xlabel('time (s)')
subplot(3,3,3)
xlabel('frequency (Hz)')
ylabel('amplitude')
title('FFT')
subplot(3,3,6)
xlabel('frequency (Hz)')
ylabel('amplitude')
Mohan kand on 21 Mar 2023
Yes to avoid the noise I used a low pass filter with matlab bulid fuction "lowpass".
In the input of lowpass fuction I use Passband frequency of 2kHz and sample rate of 2000Hz. Is it right ? See the attached plot after applied LP filter.
Literature publised used bandpass filters.

Star Strider on 21 Mar 2023
Edited: Star Strider on 21 Mar 2023
Perhaps something like this —
Fs = 2000; % Sampling Frequency (Hz)
Fn = Fs/2; % Nyquist Frequency
Ts = 1/Fs; % Sampling Time (sec), pulse repitation time
t = linspace(0, L, L)*Ts; % Time Vector (sec)
figure(1)
grid
xlabel('Time')
ylabel('Amplitude') % Signal centralization
NFFT = 2^nextpow2(L)
NFFT = 131072
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTs(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
xlim([0 20]) [envu,envl] = envelope(abs(FTs(Iv))*2, 10, 'peaks'); % Signal Envelope
[pks,locs] = findpeaks(envu, 'NPeaks',1); % Return Single Peak
cidx = find(diff(sign(envu - pks/2))); % Half-Magnitude Approximate Incides
[~,cidxx] = mink(abs(cidx-locs),2); % Half-Magnitude Approximate Indices
cidx = cidx(cidxx); % Choose Two Closest Indices To Central Peak Index
for k = 1:numel(cidx)
idxrng = max(1,cidx(k)-1) : min(numel(envu),cidx(k)+1);
xv(k) = interp1(envu(idxrng),Fv(idxrng),5E-3); % 'Precise' Frequency Values
end
figure
plot(Fv, abs(FTs(Iv))*2)
hold on
plot(Fv, envu, '-r')
hold off
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
xlim([0 10])
xline(xv,'-k', "Passband Cutoff "+xv+" Hz") % Use Frequency Values To Design Bandpass Filter % filter
Passband_Frequency_kHz = xv;
Filtered = bandpass(radarSig,Passband_Frequency_kHz,Fs, 'ImpulseResponse','iir'); % Filter Signal
[fenvu,fenvl] = envelope(Filtered, 10, 'peak'); % Envelope Of Foltered Signal
[pks,locs] = findpeaks(fenvu, 'MinPeakProminence',0.05); % Detect Upper Envelope Peaks
RRate = 1/mean(diff(t(locs)))*60; % Calculate Respiratory Rate
% FFT
% [radarSpectrumForResp_filtered, ~] = FFT(Filtered, Fs) ;
% Estimate respiratory rate
figure(1)
subplot(3,3,[1 2])
xlabel('time (s)')
ylabel('amplitude (V)')
grid
subplot(3,3,[4 5])
plot(t, Filtered)
hold on
plot(t, fenvu, 'r')
plot(t(locs), pks, 'vg', 'MarkerFaceColor','g') % Plot Peaks
hold off
ylabel('amplitude (V)')
xlabel('time (s)')
grid
text(median(t), min(ylim), sprintf('Resporatory Rate = %.2f / minute',RRate), 'Horiz','center', 'Vert','bottom') pkst = 1×9
0.8670 4.5405 12.0416 20.7287 30.6128 35.7703 42.9574 50.6129 58.1720
% subplot(3,3,3)
% xlabel('frequency (Hz)')
% ylabel('amplitude')
% title('FFT')
% subplot(3,3,6)
% xlabel('frequency (Hz)')
% ylabel('amplitude')
The envelope of the filtered signal may provide additional useful information.
EDIT — (21 Mar 2023 at 20:44)
Added code to detect upper envelope peaks and calculate respiratory rate.
.
Star Strider on 22 Mar 2023
As always, my pleasure!
Would it be possible to ask you something regarding your code if I don't understand anything? Thanks again !
Yes!
.

Mohan kand on 23 Mar 2023
Edited: Mohan kand on 23 Mar 2023
@Star Strider Can you explain how can we decide 'MinPeakProminence' value. When I use a value of 0.01 the RR is 23.01 and when I use a value of 0.05 than the RR is 21.91 (see the attached figures). Which one will be the correct value ? I check matlab help regarding this but cannot understand clearly.
Also please explain line PksThrshld = pks*0.3, why you use this ?
Thanks !
Star Strider on 27 Mar 2023
It would appear to me to be acceptable, however the filter appears to me to be a bit ‘aggressive’ and eliminates much of the signal detail. I made some changes (subtracting the mean of the signals before using fft on them) and co-plotted the time domain signals before and after filltering.
Fs = 2000; % Sampling Frequency (Hz)
Fn = Fs/2; % Nyquist Frequency
Ts = 1/Fs; % Sampling Time (sec), pulse repitation time
t = linspace(0, L, L)*Ts; % Time Vector (sec)
% fft of raw data
NFFT = 2^nextpow2(L)
NFFT = 65536
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure(2)
subplot (2,2,1)
plot(Fv, abs(FTs(Iv))*2)
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('FFT signal original')
% fft of filtered data
FTs_filtered = fft((FF_sg-mean(FF_sg)).*hann(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
subplot (2,2,2)
plot(Fv, abs(FTs_filtered(Iv))*2)
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('FFT signal SG_filtered')
% fft of the filtered data by the fft of the raw data
tt = abs(FTs_filtered./FTs);
subplot (2,2,3)
plot(Fv, tt(Iv)*2)
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('FFT signal after divide') figure
hold on
plot(t, FF_sg, 'DisplayName','SG Filtered')
hold off
grid
xlim([min(t) max(t)])
xlabel('Time (s)')
ylabel('SG Filtered Signal') I suggest experimenting with a shorter value for ‘framelen’ since it is not obvious to me how to get RR information from the filtered signal as it currently exists. Too much detail is lost, in my opinion.
.

Mohan kand on 27 Mar 2023
@Star Strider really appreciate your help ! I try to calculate the RRate based on 'tt' but getting NaN values. Could you please chek the below code and help :
[fenvu,fenvl] = envelope(tt, 10, 'peak'); % Envelope Of Foltered Signal
[pks,locs] = findpeaks(fenvu,'MinPeakProminence',max(fenvu)*3.8);% Detect Upper Envelope Peaks
% [pks,locs] = findpeaks(fenvu, 'MinPeakProminence',0.05);
RRate = 1/mean(diff(t(locs)))*60% Calculate Respiratory Rate
Star Strider on 29 Mar 2023
Do you mean on-half of "fenvu" or the filterd singnal "FF_sg" or the orginal data set?
I am referring to half of the maximum of whatever I am using as the first argument to findpeaks generally. In this instance, it is ‘fenvu’.

Mohan kand on 18 Apr 2023
##### 2 CommentsShow 1 older commentHide 1 older comment
Mohan kand on 19 Apr 2023
Edited: Mohan kand on 19 Apr 2023
@Star Strider I have uploaded a video and a signal file. I want to synchronize waveform data with an associated video. Thanks !