Main Content

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

Figure contains 6 axes objects. Axes object 1 with title Sample GDN0006Apnea7 contains an object of type line. Axes object 2 contains an object of type line. Axes object 3 with title Sample GDN0006Resting58 contains an object of type line. Axes object 4 contains an object of type line. Axes object 5 with title Sample GDN0006Valsalva113 contains an object of type line. Axes object 6 contains an object of type line.

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

Figure contains 9 axes objects. Axes object 1 with title Sample GDN0006Apnea9 contains an object of type line. Axes object 2 contains an object of type line. Axes object 3 contains an object of type line. Axes object 4 with title Sample GDN0006Resting15 contains an object of type line. Axes object 5 contains an object of type line. Axes object 6 contains an object of type line. Axes object 7 with title Sample GDN0006Valsalva62 contains an object of type line. Axes object 8 contains an object of type line. Axes object 9 contains an object of type line.

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

Figure contains 8 axes objects. Axes object 1 with title ECG Signal and Radar Signal MODWTMRA of Sample GDN0001Resting1 contains an object of type line. Axes object 2 contains an object of type line. Axes object 3 contains an object of type line. Axes object 4 contains an object of type line. Axes object 5 contains an object of type line. Axes object 6 contains an object of type line. Axes object 7 contains an object of type line. Axes object 8 contains an object of type line.

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

Figure contains an axes object. The axes object with title Training Information contains 4 objects of type line, scatter. These objects represent Training Loss of ConvAE + LSTM, Validation Loss of ConvAE + LSTM, Training Loss of ConvAE + LSTM + modwtLayer, Validation Loss of ConvAE + LSTM + modwtLayer.

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)

Figure contains 12 axes objects. Axes object 1 with title Sample GDN0006Apnea9 contains an object of type line. Axes object 2 contains an object of type line. Axes object 3 contains an object of type line. Axes object 4 contains an object of type line. Axes object 5 with title Sample GDN0006Resting15 contains an object of type line. Axes object 6 contains an object of type line. Axes object 7 contains an object of type line. Axes object 8 contains an object of type line. Axes object 9 with title Sample GDN0006Valsalva62 contains an object of type line. Axes object 10 contains an object of type line. Axes object 11 contains an object of type line. Axes object 12 contains an object of type line.

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

Figure contains an axes object. The axes object with title Distribution of Test MSEs contains 2 objects of type histogram. These objects represent 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