How to calculate respiratory rate from a Doppler radar signal

14 views (last 30 days)
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);
% Get Header
header = Data(cellfun('isempty',numericalArray));
data_ADC = vertcat(numericalArray{:});
figure(5);plot(data_ADC);
Fs = 2000; % Sampling Frequency (Hz)
Fn = Fs/2; % Nyquist Frequency
Ts = 1/Fs; % Sampling Time (sec), pulse repitation time
L = numel(data_ADC);
t = linspace(0, L, L)*Ts; % Time Vector (sec)
figure(1)
plot(t, data_ADC)
grid
xlabel('Time')
ylabel('Amplitude')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Signal centralization
radarSig = data_ADC - mean(data_ADC) ;
% filter
Passband_Frequency_kHz = 2;
Filtered = lowpass(radarSig,Passband_Frequency_kHz,Fs);
% FFT
[radarSpectrum_orignal, xf] = FFT(radarSig, Fs);
[radarSpectrumForResp_filtered, ~] = FFT(Filtered, Fs) ;
% Estimate respiratory rate
[radarRR] = estRR(radarSpectrumForResp_filtered, xf) ;
figure(1)
subplot(3,3,[1 2])
plot(t, radarSig)
xlabel('time (s)')
ylabel('amplitude (V)')
title('Radar raw signal')
subplot(3,3,[4 5])
plot(t, Filtered)
ylabel('amplitude (V)')
xlabel('time (s)')
subplot(3,3,3)
plot(xf, radarSpectrum_orignal)
xlabel('frequency (Hz)')
ylabel('amplitude')
title('FFT')
subplot(3,3,6)
plot(xf, radarSpectrumForResp_filtered)
xlabel('frequency (Hz)')
ylabel('amplitude')
title(['FFT', 'RR(radar): ', num2str(radarRR),' bpm'])
  6 Comments
Image Analyst
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
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.

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 21 Mar 2023
Edited: Star Strider on 21 Mar 2023
Perhaps something like this —
LD = load(websave('radar_signal','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1331820/radar_signal.mat'));
data_ADC = LD.data_ADC;
% figure(5);plot(data_ADC);
Fs = 2000; % Sampling Frequency (Hz)
Fn = Fs/2; % Nyquist Frequency
Ts = 1/Fs; % Sampling Time (sec), pulse repitation time
L = numel(data_ADC);
t = linspace(0, L, L)*Ts; % Time Vector (sec)
figure(1)
plot(t, data_ADC)
grid
xlabel('Time')
ylabel('Amplitude')
% Signal centralization
radarSig = data_ADC - mean(data_ADC) ;
NFFT = 2^nextpow2(L)
NFFT = 131072
FTs = fft(radarSig.*hann(L), NFFT)/L;
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
% [radarSpectrum_orignal, xf] = FFT(radarSig, Fs);
% [radarSpectrumForResp_filtered, ~] = FFT(Filtered, Fs) ;
% Estimate respiratory rate
% [radarRR] = estRR(radarSpectrumForResp_filtered, xf) ;
figure(1)
subplot(3,3,[1 2])
plot(t, radarSig)
xlabel('time (s)')
ylabel('amplitude (V)')
title('Radar raw signal')
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)
% plot(xf, radarSpectrum_orignal)
% xlabel('frequency (Hz)')
% ylabel('amplitude')
% title('FFT')
% subplot(3,3,6)
% plot(xf, radarSpectrumForResp_filtered)
% xlabel('frequency (Hz)')
% ylabel('amplitude')
% title(['FFT', 'RR(radar): ', num2str(radarRR),' bpm'])
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.
.
  6 Comments
Mohan kand
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
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!
.

Sign in to comment.

More Answers (3)

Mohan kand
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 !
  7 Comments
Mohan kand
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)
FTs = fft(radarSig.*hann(L), NFFT)/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
FF_sg = sgolayfilt(radarSig,3,71);
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
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.
LD = load(websave('data_for SG','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1337354/data_for%20SG.mat'));
data_ADC = LD.data_ADC;
radarSig = data_ADC;
% figure(5);plot(data_ADC);
Fs = 2000; % Sampling Frequency (Hz)
Fn = Fs/2; % Nyquist Frequency
Ts = 1/Fs; % Sampling Time (sec), pulse repitation time
L = numel(data_ADC);
t = linspace(0, L, L)*Ts; % Time Vector (sec)
% fft of raw data
NFFT = 2^nextpow2(L)
NFFT = 65536
FTs = fft((radarSig-mean(radarSig)).*hann(L), NFFT)/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')
xlim([0 10]) % ADDED
% fft of filtered data
FF_sg = sgolayfilt(radarSig,3,71);
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')
xlim([0 10]) % ADDED
% 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')
xlim([0 10]) % ADDED
figure
plot(t, radarSig, 'DisplayName','Original')
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.
.

Sign in to comment.


Mohan kand
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
  13 Comments
Mohan kand
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
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’.

Sign in to comment.


Mohan kand
Mohan kand on 18 Apr 2023
  2 Comments
Star Strider
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
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 !
https://ch.mathworks.com/matlabcentral/answers/1949023-how-can-i-synchronize-radar-signal-with-video-signal?s_tid=prof_contriblnk

Sign in to comment.

Categories

Find more on Signal Generation and Preprocessing in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!