# How to calculate respiratory rate from a Doppler radar signal

23 views (last 30 days)

Show older comments

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
on 21 Mar 2023

### Accepted Answer

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)

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')

% 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

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!

.

### More Answers (3)

Mohan kand
on 23 Mar 2023

Edited: Mohan kand
on 23 Mar 2023

##### 7 Comments

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)

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.

.

Mohan kand
on 27 Mar 2023

##### 13 Comments

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 Comments

Star Strider
on 18 Apr 2023

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!