How to calculate respiratory rate from a Doppler radar signal

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

[Image Analyst]

A screenshot of your signal plotted would help.

[Mohan kand]

Please see the screenshot of the signal

[Image Analyst]

Did you try using sgolayfilt followed by findpeaks?

[Mohan kand]

No I did't use sgolayfit. Actually I cannot decide what kind of filter I need to use (low pass filter or high pass filter) before applying FFT. In the matlab code I use lowpass filter.

[Image Analyst]

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]

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]

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.

.

[Mohan kand]

@Star Strider Thanks for the detail code, It is really helpful. Now I can use the Filtered signal to calculate the Respiratory rate ?

[Star Strider]

As always, my pleasure!

I added these lines:

[pks,locs] = findpeaks(fenvu, 'MinPeakProminence',0.05); % Detect Upper Envelope Peaks

RRate = 1/mean(diff(t(locs)))*60; % Calculate Respiratory Rate

and a text call to display it on the plot.

That should work. The 'MinPeakProminence' value may need to be tweaked slightly for other data sets. .

.

[Mohan kand]

@Star Strider when I using the code for other data (see attached mat file) the xv matrix producing [NaN NaN]. Do I need to change something in the for loop ?

[Star Strider]

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]

@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]

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]

@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]

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

.

[Mohan kand]

@Star Strider when the filtered signal showes similar pattarn link raw signal than the result of Respiratory rate is fine but when the filtered signal pattern is different from the raw signal than it is showing incorrect results. See the attached files. Same thing is happning with the Data Set-2 (you used) but Data Set-1 is fine.

Any idea need to fine tune something? thanks !

Correct RR:

Incorrect RR:

[Star Strider]

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]

@Star Strider Actually I used Google fit app to detect the RR and compare with our matlab code and noticed that difference. I don't know how accurate is Google fit to monitor RR. Anyway I will check more about it.

Regarding 'sgolayfilt' you mean to apply sgolayfilt on the rawdata or FFT data (I our case radarSig or Fv) and than use findpeaks ? Thanks !

[Star Strider]

I doubt that I would trust the Google fit app. If you do not have an independent measure of respiratory rate (there are several ways to do this, direct spirometry being the most reliable), it may be difficult to determine what the actual radar respiratory signal is and respiratory rate are.

Google has absolutely no idea! It’s guessing just as I am, without my equivalent medical knowledge.

Apply sgolayfilt on the raw ‘radarSig’ data. You can compare the fft of the filtered data to the fft of the raw data to see the filtering effect. If you divide element-wise the fft of the filtered data by the fft of the raw data you will get the effective transfer function of the sgolayfilt filter. That will help you understand how it works.

.

[Mohan kand]

@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]

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]

@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]

The ‘tt’ vector in your latest code is the m,agnitude value of the Fourier transform of the transfer function of the Savitzky-Golay filter output.

You probably want to use ‘FF_sg’ instead.

[Mohan kand]

Still getting NaN value for RRate. Do I have to change the MinPeakProminence ?

[fenvu,fenvl] = envelope(FF_sg, 10, 'peak'); % Envelope Of Foltered Signal

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

RRate = 1/mean(diff(t(locs)))*60% Calculate Respiratory Rate

[Star Strider]

Please post the entire code. I have no idea where that belongs in the latest version.

[Mohan kand]

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]

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.

.

[Mohan kand]

Thank you so much but it really helped me a lot. I will use reference method to compare the results. Also you are right the rate is too fast because I tested the sensor on myself with fast breathing.

[Star Strider]

As always, my pleasure!

[Mohan kand]

@Star Strider could you please explain why you devide by 3.8 ? This value will be fix for all data ?

[pks,locs] = findpeaks(fenvu,'MinPeakProminence',max(fenvu)/3.8);

Also please expain how can I determine the window length in envelope?

[fenvu,fenvl] = envelope(FF_sg, 650, 'peak'); % Envelope Of Foltered Signal

[Star Strider]

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.

.

[Mohan kand]

Also please expain how can I determine the window length in envelope? [fenvu,fenvl] = envelope(FF_sg, 650, 'peak'); % Envelope Of Foltered Signal

[Star Strider]

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]

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]

‘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]

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 !

[Star Strider]

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]

@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