# How to calculate respiratory rate from a Doppler radar signal

23 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')
Image Analyst on 21 Mar 2023
A high pass filter would just make it noisier. You'd want a low pass filter. A Savitzky-Golay filter does a sliding window fit of a polynomial to your signal, which is essentially a non-linear low pass filter in the time domain. Doing a linear low pass filter in the Fourier/spectral domain by zeroing out higher temporal frequencies and then inverse transforming would give a smoother looking signal. It might be just a case of try it and see. Try both ways while adjusting parameters. What kind of filter does the literature (published papers) say to use?
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.
.
Mohan kand on 22 Mar 2023
@Star Strider Really appriciate for your help. Now It is working for all data type. I will play with MinPeakProminence value to calculate the RR.
Would it be possible to ask you something regarding your code if I don't understand anything? Thanks again !
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 !
Mohan kand on 27 Mar 2023
@Star Strider as you suggested I have applied following code, Please check if I am doing right, I have applied 'frameleng of 71'. I have added the matfile.
% fft of raw data
NFFT = 2^nextpow2(L)
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.*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)./abs(FTs);
subplot (2,2,3)
plot(Fv, tt(Iv)*2)
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('FFT signal after divide')
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
Mohan kand on 29 Mar 2023
As you said: "I usually start with defining the ‘MinPeakProminence’ value as one-half the maximum in the data set, and then experiment". Do you mean on-half of "fenvu" or the filterd singnal "FF_sg" or the orginal data set?
[fenvu,fenvl] = envelope(FF_sg);
[pks,locs] = findpeaks(fenvu,'MinPeakProminence',max(fenvu)/3.8);
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
Star Strider on 18 Apr 2023
I would at least need the signals in order to attempt that. I have no idea what the video is, and I am not certain there would be a way to get data from it to allow the synchronization.
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 !