How to calculate respiratory rate from a Doppler radar signal

Question

Mohan kand 2023-3-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1932650-how-to-calculate-respiratory-rate-from-a-doppler-radar-signal

编辑： Mohan kand 2023-4-19

采纳的回答： Star Strider

在 MATLAB Online 中打开

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 个评论
显示 4更早的评论隐藏 4更早的评论

Image Analyst 2023-3-21

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 2023-3-21

LP_filter.PNG

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.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-3-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1932650-how-to-calculate-respiratory-rate-from-a-doppler-radar-signal#answer_1197990

编辑：Star Strider 2023-3-21

在 MATLAB Online 中打开

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 个评论
显示 4更早的评论隐藏 4更早的评论

Star Strider 2023-3-22

在 MATLAB Online 中打开

I needed to make my code a bit more robust to handle more conditions in the data. I tested this with both data sets (using the ‘k1’ loop), and it appears to work well with both.

Try this —

LD{1} = load(websave('radar_signal','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1331820/radar_signal.mat'));

LD{2} = load(websave('data_new','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1332695/data_new.mat'));

for k1 = 1:numel(LD)

fprintf(['\n' repmat('\t',1,5) '— Data Set %d —\n'],k1)

data_ADC = LD{k1}.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

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,wdth,prm] = findpeaks(envu, 'NPeaks',1); % Return Single Peak

pks = pks(1); locs = locs(1); wdth = wdth(1); prm = prm(1);

PksThrshld = pks*0.3; % Threshold For Peak Width Coordinates (Defines 'bandpass' Filter Passband)

cidx = find(diff(sign(envu - PksThrshld))); % Half-Magnitude Approximate Incides

cidxx = find(cidx < locs); % Half-Magnitude Approximate Indices

cidx = cidx(max(cidxx)+[0 1]); % 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),PksThrshld); % 'Precise' Frequency Values

end

figure

plot(Fv, abs(FTs(Iv))*2)

hold on

plot(Fv, envu, '-r')

plot(Fv(locs), pks, 'pg')

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

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

sgtitle("Data Set "+k1)

end

— Data Set 1 —

NFFT = 131072

— Data Set 2 —

NFFT = 131072

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

This should be reasonably robust, however it may need tweaking for data sets that are significaantly different from these.

.

Mohan kand 2023-3-22

@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 2023-3-22

As always, my pleasure!

‘Would it be possible to ask you something regarding your code if I don't understand anything? Thanks again !’

Yes!

.

请先登录，再进行评论。

Answer 2

Mohan kand 2023-3-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1932650-how-to-calculate-respiratory-rate-from-a-doppler-radar-signal#answer_1199170

编辑：Mohan kand 2023-3-23

@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 个评论
显示 5更早的评论隐藏 5更早的评论

Star Strider 2023-3-23

在 MATLAB Online 中打开

The 'MinPeakProminence' value can be difficult to define for any specific data set, however other approaches, such as 'MinPeakHeignt' and others have even more problems because of the nature of the data. I decided just now to define it as:

[pks,locs,wdth,prm] = findpeaks(fenvu, 'MinPeakProminence',max(envu)*3.8); % Detect Upper Envelope Peaks

defining it as a function of the maximum heak height. That seems to work rreasonably well on these two data sets. It is difficult to define precisely, even relative to the individual data sets, because of the varying nature of the data.

The ‘PksThreshld’ is defined in its comment:

PksThrshld = pks*0.3; % Threshold For Peak Width Coordinates (Defines 'bandpass' Filter Passband)

It is the value of the main peak relative to the main peak height of the fft result that I choose to define the passbands of the bandpass filter, so the frequencies on either side of the main peak that correspond to that value. Those frequencies are later determined precisely as ‘xv’. The passbands are chosen to provide as much detail as necessasry in the filtered signal, and to filter out any baselilne variation and high-frequency noise.

The tweaked code —

LD{1} = load(websave('radar_signal','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1331820/radar_signal.mat'));

LD{2} = load(websave('data_new','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1332695/data_new.mat'));

for k1 = 1:numel(LD)

fprintf(['\n' repmat('\t',1,5) '— Data Set %d —\n'],k1)

data_ADC = LD{k1}.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

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,wdth,prm] = findpeaks(envu, 'NPeaks',1); % Return Single Peak

pks = pks(1); locs = locs(1); wdth = wdth(1); prm = prm(1);

PksThrshld = pks*0.3; % Threshold For Peak Width Coordinates (Defines 'bandpass' Filter Passband)

cidx = find(diff(sign(envu - PksThrshld))); % Half-Magnitude Approximate Incides

cidxx = find(cidx < locs); % Half-Magnitude Approximate Indices

cidx = cidx(max(cidxx)+[0 1]); % 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),PksThrshld); % 'Precise' Frequency Values

end

figure

plot(Fv, abs(FTs(Iv))*2)

hold on

plot(Fv, envu, '-r')

plot(Fv(locs), pks, 'pg')

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,wdth,prm] = findpeaks(fenvu, 'MinPeakProminence',max(envu)*3.8); % 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

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

sgtitle("Data Set "+k1)

end

— Data Set 1 —

NFFT = 131072

— Data Set 2 —

NFFT = 131072

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

.

Star Strider 2023-3-23

I have no idea what the ‘actual’ respiratory rate is in any given data set. A resting respiratory rate of 15 would seem normal to me, a resporatory rate of 30 would indicate tachypnea, although I would have to know the context. Experiment with the various parameter settings (specifically MinPeakProminence' that would seem to govern how the peaks are defined) to get the appropriate result. There is nothing in the (rather noisy) data that would point to a specific value for that. I was using the wwaveforms in the reference ass a guide as to what the results should look like.

Another option to consider instead of the frequency-selective filter that I use here, is to use a Savitzky-Golay filter (sgolayfilt) to do the filtering and return the appropriate result. (Use it in place of the frequency-selective filter I designed here.) I usually use an order 3 polynomial, and then ajdust the ‘framelen’ value to get the result I want. It essenmtially implements a FIR filter on the data.

I went back and looked at the reference. They do not include their MATLAB code (at least that I could see), so I have no idea what signal processing design choices they made to get their results.

.

Mohan kand 2023-3-27

data_for SG.mat

@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 2023-3-27

在 MATLAB Online 中打开

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.

.

请先登录，再进行评论。

Answer 3

Mohan kand 2023-3-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1932650-how-to-calculate-respiratory-rate-from-a-doppler-radar-signal#answer_1201889

@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 个评论
显示 11更早的评论隐藏 11更早的评论

Mohan kand 2023-3-27

Here is the entire code

clear all;

close all;

[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{:});

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)

% Signal centralization

radarSig = data_ADC - mean(data_ADC);

% 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])

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

[fenvu,fenvl] = envelope(FF_sg, 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 2023-3-27

在 MATLAB Online 中打开

There are still problems, not the least of which is that we do not know what the actual respirations are because there are no actual spirometric data to compare with the radar signal to standardise it.

% % 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));

% % radarSig = data_ADC;

% % data_ADC = vertcat(numericalArray{:});

LD = load(websave('data_for SG','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1337354/data_for%20SG.mat'));

data_ADC = LD.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)

% Signal centralization

radarSig = data_ADC - mean(data_ADC);

% 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])

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

[fenvu,fenvl] = envelope(FF_sg, 650, '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

RRate = 35.7075

figure

% plot(t, radarSig, 'DisplayName','Original')

hold on

plot(t, FF_sg, 'DisplayName','SG Filtered')

plot(t, fenvu, '-r')

plot(t(locs), pks, 'vr', 'MarkerFaceColor','r')

hold off

grid

xlim([min(t) max(t)])

xlabel('Time (s)')

ylabel('SG Filtered Signal')

Here, I changed the envelope window argument and the findpeaks 'MinPeakProminence' value. You can refine these further, as well as the sgolayfilt ‘framelen’ values to get a reasonably representative respiratory signal. I can tell by inspection that the current analysis is not correct in that it neither produces the correct waveforms nor an appropriate result. The rate is much too fast, and the rhythm is too irregular for a healthy resting signal (although it may be appropriate in the context of a sick patient, however I have no idea what that context is).

Again, you really need some sort of objective measure, such as a spirometric signal, to use as a standard reference. Without it, getting the correct result simply amounts to guessing.

.

Star Strider 2023-3-28

Actually, not specifically. The ‘3.8’ value just happened to be there, so that value can be anything, really, and if that value had not worked, I would have changed it.

I usually start with defining the ‘MinPeakProminence’ value as one-half the maximum in the data set, and then experiment. Another option is to return all the prominence values and then use the maxk function to find the 5 or so largest, and then use those to estimate the correct value. I generally need to define a value that will adapt to different data sets in the same series (such as yours), and taking a fraction of the maximum is usually a good initial approach. I like using 'MinPeakProminence' because it is the most adaptable of all the findpeaks arguments.

I wish that it was possible to figure out a way to have findpeaks and other necessary functions automatically adapt to each data set and always do exactly what I want them to do. If something like that exists, I have not yet discovered it, although I am always looking for it.

.

Star Strider 2023-3-28

A lot of digital signal processing is entirely heruistic — has to be manually adapted to a particular situation with the best result decided by trial and error. Sometimes, this can be helped significantly by analysing the signal with Fouirier transforms or other techniques, as I did when creating the bandpass filter. In this instance, since it is difficult to determine the best choice with these data, I experiment until I get the sort of result I want. Even with the same data, that might change between the unfiltered original data and filtered data, since filtered data would likely have a different frequency content than the original data, and so have a different window for envelope. If the sampling frrequency is different, that also can complicate the parameter value choice (although the sampling frequency seems to be constant with these data).

So experiment until it looks like the result you want. That parameter choice will also likely work with other signals if they have the same sampling frequency and if you do the same filtering and other pre-processing on them.

Mohan kand 2023-3-29

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 2023-3-29

‘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’.

请先登录，再进行评论。

Answer 4

Mohan kand 2023-4-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1932650-how-to-calculate-respiratory-rate-from-a-doppler-radar-signal#answer_1218143

hi @Star Strider any idea about my question below :

https://ch.mathworks.com/matlabcentral/answers/1949023-how-can-i-synchronize-radar-signal-with-video-signal?s_tid=prof_contriblnk

Thanks !

2 个评论
显示无隐藏无

Star Strider 2023-4-18

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 2023-4-19

编辑：Mohan kand 2023-4-19

@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

请先登录，再进行评论。

How to calculate respiratory rate from a Doppler radar signal

6 个评论
显示 4更早的评论隐藏 4更早的评论

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（3 个）

7 个评论
显示 5更早的评论隐藏 5更早的评论

13 个评论
显示 11更早的评论隐藏 11更早的评论

2 个评论
显示无隐藏无

另请参阅

类别

标签

Community Treasure Hunt

How to calculate respiratory rate from a Doppler radar signal

6 个评论 显示 4更早的评论隐藏 4更早的评论

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（3 个）

7 个评论 显示 5更早的评论隐藏 5更早的评论

13 个评论 显示 11更早的评论隐藏 11更早的评论

2 个评论 显示 无隐藏 无

另请参阅

类别

标签

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论

13 个评论
显示 11更早的评论隐藏 11更早的评论

2 个评论
显示无隐藏无