Human Health Monitoring Using Continuous Wave Radar and Deep Learning
This example shows how to reconstruct electrocardiogram (ECG) signals via continuous wave (CW) radar signals using deep learning neural networks.
Radar is now being used for vital sign monitoring. This method offers many advantages over wearable devices. It allows non-contact measurement which is preferred for use in cases of daily use of long-term monitoring. However, the challenge is that we need to convert radar signals to vital signs or to meaningful biosignals that can be interpreted by physicians. Current traditional methods based on signal transform and correlation analysis can capture periodic heartbeats but fail to reconstruct the ECG signal from the radar returns. This example shows how AI, specifically a deep learning network, can be used to reconstruct ECG signals solely from the radar measurements.
This example uses a hybrid convolutional autoencoder and bidirectional long short-term memory (BiLSTM) network as the model. Then, a wavelet multiresolution decomposition layer, maximal overlap discrete wavelet transform (MODWT) layer, is introduced to improve the performance. The example compares the network using a 1-D convolutional layer and network using a MODWT layer.
Data Description
The dataset [1] presented in this example consists of synchronized data from a CW radar and ECG signals measured simultaneously by a reference device on 30 healthy subjects. The implemented CW radar system is based on the Six-Port technology and operates at 24 GHz in the Industrial Scientific and Medical (ISM) band.
Due to the large volume of the original dataset, for efficiency of model training, only a small subset of the data is used in this example. Specifically, the data from three scenarios, resting, apnea, and Valsalva maneuver, is selected. Further, the data from subjects 1-5 is used to train and validate the model. The data from subject 6 is used to test the trained model.
Also, because the main information contained in the ECG signal is usually located in a frequency band less than 100 Hz, all signals are downsampled to 200 Hz and divided into segments of 1024 points, i.e. signals of approximately 5s.
Download and Prepare Data
The data has been uploaded to this location: https://ssd.mathworks.com/supportfiles/SPT/data/SynchronizedRadarECGData.zip.
Download the dataset using the downloadSupportFile function. The whole dataset is approximately 16 MB in size. It contains two folders, trainVal for training and validation data and test for test data. Inside each of them, ECG signals and radar signals are stored in two separate folders, ecg and radar.
datasetZipFile = matlab.internal.examples.downloadSupportFile('SPT','data/SynchronizedRadarECGData.zip'); datasetFolder = fullfile(fileparts(datasetZipFile),'SynchronizedRadarECGData'); if ~exist(datasetFolder,'dir') unzip(datasetZipFile,datasetFolder); end
Create signal datastores to access the data in the files.
radarTrainValDs = signalDatastore(fullfile(datasetFolder,"trainVal","radar")); radarTestDs = signalDatastore(fullfile(datasetFolder,"test","radar")); ecgTrainValDs = signalDatastore(fullfile(datasetFolder,"trainVal","ecg")); ecgTestDs = signalDatastore(fullfile(datasetFolder,"test","ecg"));
View the categories and distribution of the data contained in the training and test sets. Note the GDN000X represents measurement data from subject X and not every subject has data for all three scenarios.
trainCats = filenames2labels(radarTrainValDs,'ExtractBefore','_radar'); summary(trainCats)
GDN0001_Resting 59
GDN0001_Valsalva 97
GDN0002_Resting 60
GDN0002_Valsalva 97
GDN0003_Resting 58
GDN0003_Valsalva 103
GDN0004_Apnea 14
GDN0004_Resting 58
GDN0004_Valsalva 106
GDN0005_Apnea 14
GDN0005_Resting 59
GDN0005_Valsalva 105
testCats = filenames2labels(radarTestDs,'ExtractBefore','_radar'); summary(testCats)
GDN0006_Apnea 14
GDN0006_Resting 59
GDN0006_Valsalva 127
Apply normalization on ECG signals. Center each signal by subtracting its median and rescale it so that its maximum peak is 1.
ecgTrainValDs = transform(ecgTrainValDs,@helperNormalize); ecgTestDs = transform(ecgTestDs,@helperNormalize);
Combine radar and ECG signal datastores. Then, further split the training and validation dataset. Use 90% of data for training and 10% for validation. Set the random seed so that data segmentation and visualization results are reproducible.
trainValDs = combine(radarTrainValDs,ecgTrainValDs);
testDs = combine(radarTestDs,ecgTestDs);
rng("default")
splitIndices = splitlabels(trainCats,0.90);
numTrain = length(splitIndices{1})numTrain = 747
numVal = length(splitIndices{2})numVal = 83
numTest = length(testCats)
numTest = 200
trainDs = subset(trainValDs,splitIndices{1});
valDs = subset(trainValDs,splitIndices{2});Because the dataset used here is not large, read the training, testing, and validation data into memory.
trainData = readall(trainDs); valData = readall(valDs); testData = readall(testDs);
Preview Data
Plot a representative of each type of signal. Notice that it is almost impossible to identify any correlation between the radar signals and the corresponding reference ECG measurements.
numCats = cumsum(countcats(testCats)); previewindices = [randi([1,numCats(1)]),randi([numCats(1)+1,numCats(2)]),randi([numCats(2)+1,numCats(3)])]; helperPlotData(testDs,previewindices);
Train Convolutional Autoencoder and BiLSTM Model
Build a hybrid convolutional autoencoder and BiLSTM network to reconstruct ECG signals. The first 1-D convolutional layer filters the signal. Then, the convolutional autoencoder eliminates most of the high-frequency noise and captures the high-level patterns of the whole signal. The subsequent BiLSTM layer further finely shapes the signal details.
layers1 = [
sequenceInputLayer(1,MinLength = 1024)
convolution1dLayer(4,3,Padding="same",Stride=1)
convolution1dLayer(64,8,Padding="same",Stride=8)
batchNormalizationLayer()
tanhLayer
maxPooling1dLayer(2,Padding="same")
convolution1dLayer(32,8,Padding="same",Stride=4)
batchNormalizationLayer
tanhLayer
maxPooling1dLayer(2,Padding="same")
transposedConv1dLayer(32,8,Cropping="same",Stride=4)
tanhLayer
transposedConv1dLayer(64,8,Cropping="same",Stride=8)
tanhLayer
bilstmLayer(8)
fullyConnectedLayer(8)
dropoutLayer(0.2)
fullyConnectedLayer(4)
dropoutLayer(0.2)
fullyConnectedLayer(1)
regressionLayer];Define the training option parameters: use an Adam optimizer and choose to shuffle the data at every epoch. Also, specify radarVal and ecgVal as the source for the validation data. Use the trainNetwork function to train the model. At the same time, the training information is recorded, which will be used for performance analysis and comparison later.
Train on a GPU if one is available. Using a GPU requires Parallel Computing Toolbox™ and a supported GPU device. For information on supported devices, see GPU Computing Requirements (Parallel Computing Toolbox) (Parallel Computing Toolbox).
options = trainingOptions("adam",... MaxEpochs=600,... MiniBatchSize=600,... InitialLearnRate=0.001,... ValidationData={valData(:,1),valData(:,2)},... ValidationFrequency=100,... VerboseFrequency=100,... Verbose=1, ... Shuffle="every-epoch",... Plots="none", ... DispatchInBackground=true); [net1,info1] = trainNetwork(trainData(:,1),trainData(:,2),layers1,options);
Training on single GPU. |======================================================================================================================| | Epoch | Iteration | Time Elapsed | Mini-batch | Validation | Mini-batch | Validation | Base Learning | | | | (hh:mm:ss) | RMSE | RMSE | Loss | Loss | Rate | |======================================================================================================================| | 1 | 1 | 00:00:01 | 0.19 | 0.18 | 0.0180 | 0.0171 | 0.0010 | | 100 | 100 | 00:00:57 | 0.18 | 0.18 | 0.0166 | 0.0163 | 0.0010 | | 200 | 200 | 00:01:54 | 0.18 | 0.18 | 0.0165 | 0.0164 | 0.0010 | | 300 | 300 | 00:02:50 | 0.18 | 0.18 | 0.0161 | 0.0160 | 0.0010 | | 400 | 400 | 00:03:46 | 0.17 | 0.17 | 0.0151 | 0.0148 | 0.0010 | | 500 | 500 | 00:04:42 | 0.17 | 0.17 | 0.0144 | 0.0140 | 0.0010 | | 600 | 600 | 00:05:38 | 0.17 | 0.16 | 0.0138 | 0.0133 | 0.0010 | |======================================================================================================================| Training finished: Max epochs completed.
Analyze Performance of Trained Model
Randomly pick a representative sample of each type from the test dataset to visualize and get an initial intuition about the accuracy of the reconstructions of the trained model.
testindices = [randi([1,numCats(1)]),randi([numCats(1)+1,numCats(2)]),randi([numCats(2)+1,numCats(3)])]; helperPlotData(testDs,testindices,net1);
Comparing the measured and reconstructed ECG signals, it can be seen that the model has been able to initially learn some correlations between the radar and ECG signals. But the results are not very satisfactory. Some peaks are not aligned with the actual peaks and the waveform shapes do not resemble those of the measured ECG. A few peaks are even completely lost.
Improve Performance Using Multiresolution Analysis and MODWT Layer
Feature extraction is often used to capture the key information of the signals, reduce the dimensionality and redundancy of the data, and help the model achieve better results. Considering that the effective information of the ECG signal exists in a certain frequency range. Use MODWT to decompose the radar signals and get the multiresolution analysis (MRA) of it as the feature.
It can be found that some components of the radar signal decomposed by MODWTMRA, such as the component of level 4, have similar periodic patterns with the measured ECG signal. Meanwhile, some components contain almost complete noise. Inspired by this, introducing MODWT layer into the model and selecting only a few level components may help the network focus more on correlated information, while also reducing irrelevant interference.
ds = subset(trainDs,1);
[~,name,~] = fileparts(ds.UnderlyingDatastores{1}.Files{1});
data = read(ds);
radar = data{1};
ecg = data{2};
levs = 1:6;
idx = 100:800;
m = modwt(radar,'sym2',max(levs));
nplot = length(levs)+2;
mra = modwtmra(m);
figure
tiledlayout(nplot,1)
nexttile
plot(ecg(idx))
title(["ECG Signal and Radar Signal MODWTMRA", "of Sample " + regexprep(name, {'_','radar'}, '')])
ylabel("Measured ECG")
grid on
d = 1;
for i = levs
d = d+1;
nexttile
plot(mra(i,idx))
ylabel(["Radar", "MODWTMRA", "Level " + i'])
grid on
end
nexttile
plot(mra(i+1,idx))
ylabel(["Radar", "MODWTMRA","Scaling","Coefficients"])
set(gcf,'Position',[0 0 700,800])
Replace the first convolution1dLayer with modwtLayer. The MODWT layer has been configured to have the same filter size and number of output channels to preserve the number of learning parameters. Based on the observations before, only components of a specific frequency range are preserved, i.e. level 3 to 5, which effectively removes unnecessary signal information that is irrelevant to the ECG reconstruction. Refer to modwtLayer (Wavelet Toolbox) documentation for more details on modwtLayer and these parameters.
A flattenLayer is also inserted after the modwtLayer to make the subsequent convolution1dLayer convolve along the time dimension, and to make the output format compatible with the subsequent bilstmLayer.
layers2 = [
sequenceInputLayer(1,MinLength = 1024)
modwtLayer('Level',5,'IncludeLowpass',false,'SelectedLevels',3:5,"Wavelet","sym2")
flattenLayer
convolution1dLayer(64,8,Padding="same",Stride=8)
batchNormalizationLayer()
tanhLayer
maxPooling1dLayer(2,Padding="same")
convolution1dLayer(32,8,Padding="same",Stride=4)
batchNormalizationLayer
tanhLayer
maxPooling1dLayer(2,Padding="same")
transposedConv1dLayer(32,8,Cropping="same",Stride=4)
tanhLayer
transposedConv1dLayer(64,8,Cropping="same",Stride=8)
tanhLayer
bilstmLayer(8)
fullyConnectedLayer(8)
dropoutLayer(0.2)
fullyConnectedLayer(4)
dropoutLayer(0.2)
fullyConnectedLayer(1)
regressionLayer];Use the same training options as before.
options = trainingOptions("adam",... MaxEpochs=600,... MiniBatchSize=600,... InitialLearnRate=0.001,... ValidationData={valData(:,1),valData(:,2)},... ValidationFrequency=100,... VerboseFrequency=100,... Verbose=1, ... Shuffle="every-epoch",... Plots="none", ... DispatchInBackground=true); [net2,info2] = trainNetwork(trainData(:,1),trainData(:,2),layers2,options);
Training on single GPU. |======================================================================================================================| | Epoch | Iteration | Time Elapsed | Mini-batch | Validation | Mini-batch | Validation | Base Learning | | | | (hh:mm:ss) | RMSE | RMSE | Loss | Loss | Rate | |======================================================================================================================| | 1 | 1 | 00:00:01 | 0.19 | 0.19 | 0.0180 | 0.0172 | 0.0010 | | 100 | 100 | 00:01:10 | 0.15 | 0.14 | 0.0112 | 0.0100 | 0.0010 | | 200 | 200 | 00:02:19 | 0.12 | 0.11 | 0.0074 | 0.0064 | 0.0010 | | 300 | 300 | 00:03:28 | 0.11 | 0.11 | 0.0063 | 0.0061 | 0.0010 | | 400 | 400 | 00:04:37 | 0.11 | 0.11 | 0.0058 | 0.0059 | 0.0010 | | 500 | 500 | 00:05:46 | 0.10 | 0.11 | 0.0053 | 0.0059 | 0.0010 | | 600 | 600 | 00:06:57 | 0.10 | 0.11 | 0.0051 | 0.0061 | 0.0010 | |======================================================================================================================| Training finished: Max epochs completed.
Compare the training and validation losses of two models. Both the training loss and validation loss of the model with MODWT layer drop much faster and more smoothly.
figure plot(info1.TrainingLoss) hold on scatter(1:length(info1.ValidationLoss),info1.ValidationLoss) plot(info2.TrainingLoss) scatter(1:length(info2.ValidationLoss),info2.ValidationLoss) hold off legend(["Training Loss of ConvAE + LSTM", ... "Validation Loss of ConvAE + LSTM", ... "Training Loss of ConvAE + LSTM + modwtLayer",... "Validation Loss of ConvAE + LSTM + modwtLayer"],"Location","eastoutside") xlabel("Epoch") title("Training Information") set(gcf,'Position',[0 0 1000,500]);
Further compare the reconstructed signals on the same test data samples. The model with modwtLayer can capture the peak position, magnitude, and shape of ECG signals very well in resting and valsalva scenarios. Even in apnea scenarios, with relatively few training samples, it can still capture the main peaks and get reconstructed signals.
helperPlotData(testDs,testindices,net1,net2)
Compare the distributions of the reconstructed signal errors for the two models on the full test set. It further illustrates that using MODWT layer improves the accuracy of the reconstruction.
ecgTestRe1 = predict(net1,testData(:,1)); loss1 = cellfun(@mse,ecgTestRe1,testData(:,2)); ecgTestRe2 = predict(net2,testData(:,1)); loss2 = cellfun(@mse,ecgTestRe2,testData(:,2)); figure h1 = histogram(loss1); hold on h2 = histogram(loss2); hold off h1.Normalization = 'probability'; h1.BinWidth = 0.003; h2.Normalization = 'probability'; h2.BinWidth = 0.003; ylabel("Probability") xlabel("MSE between Measured and Reconstructed ECG Signals") title("Distribution of Test MSEs") legend(["Model without modwtLayer","Model with modwtLayer"])
Conclusion
This example implements a convolutional autoencoder and BiLSTM network to reconstruct ECG signals from CW radar signals. The example analyzes the performance of the model with and without a MODWT Layer. It shows that the introduction of MODWT layer improves the quality of the reconstructed ECG signals.
Reference
[1] Schellenberger, S., Shi, K., Steigleder, T. et al. A dataset of clinically recorded radar vital signs with synchronised reference sensor signals. Sci Data 7, 291 (2020). https://doi.org/10.1038/s41597-020-00629-5
Appendix -- Helper Functions
helperNormalize - this function normalizes input signals by subtracting the median and dividing by the maximum value.
function x = helperNormalize(x) % This function is only intended to support this example. It may be changed % or removed in a future release. x = x-median(x); x = {x/max(x)}; end
helperPlotData - this function plots radar and ECG signals.
function helperPlotData(DS,Indices,net1,net2) % This function is only intended to support this example. It may be changed % or removed in a future release. arguments DS Indices net1 =[] net2 = [] end fs = 200; N = numel(Indices); M = 2; if ~isempty(net1) M = M + 1; end if ~isempty(net2) M = M + 1; end tiledlayout(M, N, 'Padding', 'none', 'TileSpacing', 'compact'); for i = 1:N idx = Indices(i); ds = subset(DS,idx); [~,name,~] = fileparts(ds.UnderlyingDatastores{1}.Files{1}); data = read(ds); radar = data{1}; ecg = data{2}; t = linspace(0,length(radar)/fs,length(radar)); nexttile(i) plot(t,radar) title(["Sample",regexprep(name, {'_','radar'}, '')]) xlabel(["Radar Signal","Time (s)"]) grid on nexttile(N+i) plot(t,ecg) xlabel(["Measured ECG Signal","Time (s)"]) ylim([-0.3,1]) grid on if ~isempty(net1) nexttile(2*N+i) y = predict(net1,radar); plot(t,y) grid on ylim([-0.3,1]) xlabel(["Reconstructed ECG Signal","Time (s)"]) end if ~isempty(net2) nexttile(3*N+i) y = predict(net2,radar); hold on plot(t,y) hold off grid on ylim([-0.3,1]) xlabel(["Reconstructed ECG Signal", "with modwtLayer","Time (s)"]) end end set(gcf,'Position',[0 0 300*N,150*M]) end
See Also
Objects
Functions
You can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)