This example shows how to classify radar returns with both machine and deep learning approaches. The machine learning approach uses wavelet scattering feature extraction coupled with a support vector machine. Additionally, two deep learning approaches are illustrated: transfer learning using SqueezeNet and a Long Short-Term Memory (LSTM) recurrent neural network. Note that the data set used in this example does not require advanced techniques but the workflow is described because the techniques can be extended to more complex problems.

Target classification is an important function in modern radar systems. This example uses machine and deep learning to classify radar echoes from a cylinder and a cone. Although this example uses the synthesized I/Q samples, the workflow is applicable to real radar returns.

The next section shows how to create synthesized data to train the learning algorithms.

The following code simulates the RCS pattern of a cylinder with a radius of 1 meter and a height of 10 meters. The operating frequency of the radar is 850 MHz.

c = 3e8; fc = 850e6; [cylrcs,az,el] = rcscylinder(1,1,10,c,fc); helperTargetRCSPatternPlot(az,el,cylrcs);

The pattern can then be applied to a backscatter radar target to simulate returns from different aspects angles.

cyltgt = phased.BackscatterRadarTarget('PropagationSpeed',c,... 'OperatingFrequency',fc,'AzimuthAngles',az,'ElevationAngles',el,'RCSPattern',cylrcs);

The following plot shows how to simulate 100 returns of the cylinder over time. It is assumed that the cylinder under goes a motion that causes small vibrations around bore sight, as a result, the aspect angle changes from one sample to the next.

rng default; N = 100; az = 2*randn(1,N); el = 2*randn(1,N); cylrtn = cyltgt(ones(1,N),[az;el]); plot(mag2db(abs(cylrtn))); xlabel('Time Index') ylabel('Target Return (dB)'); title('Target Return for Cylinder');

The return of the cone can be generated similarly. To create the training set, the above process is repeated for 5 arbitrarily selected cylinder radii. In addition, for each radius, 10 motion profiles are simulated by varying the incident angle following 10 randomly generated sinusoid curve around boresight. There are 701 samples in each motion profile, so there are 701-by-50 samples. The process is repeated for the cylinder target, which results in a 701-by-100 matrix of training data with 50 cylinder and 50 cone profiles. In the test set, we use 25 cylinder and 25 cone profiles to create a 701-by-50 training set. Because of the long computation time, the training data is precomputed and loaded below.

load('RCSClassificationReturnsTraining'); load('RCSClassificationReturnsTest');

As an example, the next plot shows the return for one of the motion profiles from each shape. The plots show how the values change over time for both the incident azimuth angles and the target returns.

subplot(2,2,1) plot(cylinderAspectAngle(1,:)) ylim([-90 90]) grid on title('Cylinder Aspect Angle vs. Time'); xlabel('Time Index'); ylabel('Aspect Angle (degrees)'); subplot(2,2,3) plot(RCSReturns.Cylinder_1); ylim([-50 50]); grid on title('Cylinder Return'); xlabel('Time Index'); ylabel('Target Return (dB)'); subplot(2,2,2) plot(coneAspectAngle(1,:)); ylim([-90 90]); grid on; title('Cone Aspect Angle vs. Time'); xlabel('Time Index'); ylabel('Aspect Angle (degrees)'); subplot(2,2,4); plot(RCSReturns.Cone_1); ylim([-50 50]); grid on; title('Cone Return'); xlabel('Time Index'); ylabel('Target Return (dB)');

In the wavelet scattering feature extractor, data is propagated through a series of wavelet transforms, nonlinearities, and averaging to produce low-variance representations of time series. Wavelet time scattering yields signal representations insensitive to shifts in the input signal without sacrificing class discriminability.

The key parameters to specify in a wavelet time scattering decomposition are the scale of the time invariant, the number of wavelet transforms, and the number of wavelets per octave in each of the wavelet filter banks. In many applications, the cascade of two filter banks is sufficient to achieve good performance. In this example, we construct a wavelet time scattering decomposition with the default filter banks: 8 wavelets per octave in the first filter bank and 1 wavelet per octave in the second filter bank. The invariance scale is set to 701 samples, the length of the data.

sf = waveletScattering('SignalLength',701,'InvarianceScale',701);

Next, we obtain the scattering transforms of both the training and test sets.

sTrain = sf.featureMatrix(RCSReturns{:,:},'transform','log'); sTest = sf.featureMatrix(RCSReturnsTest{:,:},'transform','log');

For this example, use the mean of the scattering coefficients taken along each path.

TrainFeatures = squeeze(mean(sTrain,2))'; TestFeatures = squeeze(mean(sTest,2))';

Create the labels for training and learning

TrainLabels = repelem(categorical({'Cylinder','Cone'}),[50 50])'; TestLabels = repelem(categorical({'Cylinder','Cone'}),[25 25])';

Fit a support vector machine model with a quadratic kernel to the scattering features and obtain the cross-validation accuracy.

template = templateSVM('KernelFunction', 'polynomial', ... 'PolynomialOrder', 2, ... 'KernelScale', 'auto', ... 'BoxConstraint', 1, ... 'Standardize', true); classificationSVM = fitcecoc(... TrainFeatures, ... TrainLabels, ... 'Learners', template, ... 'Coding', 'onevsone', ... 'ClassNames', categorical({'Cylinder','Cone'})); partitionedModel = crossval(classificationSVM, 'KFold', 5); [validationPredictions, validationScores] = kfoldPredict(partitionedModel); validationAccuracy = (1 - kfoldLoss(partitionedModel, 'LossFun', 'ClassifError'))*100

validationAccuracy = 100

Using the trained SVM, classify the scattering features obtained from the test set.

predLabels = predict(classificationSVM,TestFeatures); accuracy = sum(predLabels == TestLabels )/numel(TestLabels)*100

accuracy = 100

Plot the confusion matrix.

figure('Units','normalized','Position',[0.2 0.2 0.5 0.5]); ccDCNN = confusionchart(TestLabels,predLabels); ccDCNN.Title = 'Confusion Chart'; ccDCNN.ColumnSummary = 'column-normalized'; ccDCNN.RowSummary = 'row-normalized';

For more complex data sets, a deep learning workflow may improve performance.

SqueezeNet is a deep convolutional neural network (CNN) trained for images in 1,000 classes as used in the ImageNet Large Scale Visual Recognition Challenge (ILSVRC). In this example, we reuse the pre-trained SqueezeNet to classify radar returns belonging to one of two classes.

Load SqueezeNet.

snet = squeezenet; snet.Layers

ans = 68x1 Layer array with layers: 1 'data' Image Input 227x227x3 images with 'zerocenter' normalization 2 'conv1' Convolution 64 3x3x3 convolutions with stride [2 2] and padding [0 0 0 0] 3 'relu_conv1' ReLU ReLU 4 'pool1' Max Pooling 3x3 max pooling with stride [2 2] and padding [0 0 0 0] 5 'fire2-squeeze1x1' Convolution 16 1x1x64 convolutions with stride [1 1] and padding [0 0 0 0] 6 'fire2-relu_squeeze1x1' ReLU ReLU 7 'fire2-expand1x1' Convolution 64 1x1x16 convolutions with stride [1 1] and padding [0 0 0 0] 8 'fire2-relu_expand1x1' ReLU ReLU 9 'fire2-expand3x3' Convolution 64 3x3x16 convolutions with stride [1 1] and padding [1 1 1 1] 10 'fire2-relu_expand3x3' ReLU ReLU 11 'fire2-concat' Depth concatenation Depth concatenation of 2 inputs 12 'fire3-squeeze1x1' Convolution 16 1x1x128 convolutions with stride [1 1] and padding [0 0 0 0] 13 'fire3-relu_squeeze1x1' ReLU ReLU 14 'fire3-expand1x1' Convolution 64 1x1x16 convolutions with stride [1 1] and padding [0 0 0 0] 15 'fire3-relu_expand1x1' ReLU ReLU 16 'fire3-expand3x3' Convolution 64 3x3x16 convolutions with stride [1 1] and padding [1 1 1 1] 17 'fire3-relu_expand3x3' ReLU ReLU 18 'fire3-concat' Depth concatenation Depth concatenation of 2 inputs 19 'pool3' Max Pooling 3x3 max pooling with stride [2 2] and padding [0 1 0 1] 20 'fire4-squeeze1x1' Convolution 32 1x1x128 convolutions with stride [1 1] and padding [0 0 0 0] 21 'fire4-relu_squeeze1x1' ReLU ReLU 22 'fire4-expand1x1' Convolution 128 1x1x32 convolutions with stride [1 1] and padding [0 0 0 0] 23 'fire4-relu_expand1x1' ReLU ReLU 24 'fire4-expand3x3' Convolution 128 3x3x32 convolutions with stride [1 1] and padding [1 1 1 1] 25 'fire4-relu_expand3x3' ReLU ReLU 26 'fire4-concat' Depth concatenation Depth concatenation of 2 inputs 27 'fire5-squeeze1x1' Convolution 32 1x1x256 convolutions with stride [1 1] and padding [0 0 0 0] 28 'fire5-relu_squeeze1x1' ReLU ReLU 29 'fire5-expand1x1' Convolution 128 1x1x32 convolutions with stride [1 1] and padding [0 0 0 0] 30 'fire5-relu_expand1x1' ReLU ReLU 31 'fire5-expand3x3' Convolution 128 3x3x32 convolutions with stride [1 1] and padding [1 1 1 1] 32 'fire5-relu_expand3x3' ReLU ReLU 33 'fire5-concat' Depth concatenation Depth concatenation of 2 inputs 34 'pool5' Max Pooling 3x3 max pooling with stride [2 2] and padding [0 1 0 1] 35 'fire6-squeeze1x1' Convolution 48 1x1x256 convolutions with stride [1 1] and padding [0 0 0 0] 36 'fire6-relu_squeeze1x1' ReLU ReLU 37 'fire6-expand1x1' Convolution 192 1x1x48 convolutions with stride [1 1] and padding [0 0 0 0] 38 'fire6-relu_expand1x1' ReLU ReLU 39 'fire6-expand3x3' Convolution 192 3x3x48 convolutions with stride [1 1] and padding [1 1 1 1] 40 'fire6-relu_expand3x3' ReLU ReLU 41 'fire6-concat' Depth concatenation Depth concatenation of 2 inputs 42 'fire7-squeeze1x1' Convolution 48 1x1x384 convolutions with stride [1 1] and padding [0 0 0 0] 43 'fire7-relu_squeeze1x1' ReLU ReLU 44 'fire7-expand1x1' Convolution 192 1x1x48 convolutions with stride [1 1] and padding [0 0 0 0] 45 'fire7-relu_expand1x1' ReLU ReLU 46 'fire7-expand3x3' Convolution 192 3x3x48 convolutions with stride [1 1] and padding [1 1 1 1] 47 'fire7-relu_expand3x3' ReLU ReLU 48 'fire7-concat' Depth concatenation Depth concatenation of 2 inputs 49 'fire8-squeeze1x1' Convolution 64 1x1x384 convolutions with stride [1 1] and padding [0 0 0 0] 50 'fire8-relu_squeeze1x1' ReLU ReLU 51 'fire8-expand1x1' Convolution 256 1x1x64 convolutions with stride [1 1] and padding [0 0 0 0] 52 'fire8-relu_expand1x1' ReLU ReLU 53 'fire8-expand3x3' Convolution 256 3x3x64 convolutions with stride [1 1] and padding [1 1 1 1] 54 'fire8-relu_expand3x3' ReLU ReLU 55 'fire8-concat' Depth concatenation Depth concatenation of 2 inputs 56 'fire9-squeeze1x1' Convolution 64 1x1x512 convolutions with stride [1 1] and padding [0 0 0 0] 57 'fire9-relu_squeeze1x1' ReLU ReLU 58 'fire9-expand1x1' Convolution 256 1x1x64 convolutions with stride [1 1] and padding [0 0 0 0] 59 'fire9-relu_expand1x1' ReLU ReLU 60 'fire9-expand3x3' Convolution 256 3x3x64 convolutions with stride [1 1] and padding [1 1 1 1] 61 'fire9-relu_expand3x3' ReLU ReLU 62 'fire9-concat' Depth concatenation Depth concatenation of 2 inputs 63 'drop9' Dropout 50% dropout 64 'conv10' Convolution 1000 1x1x512 convolutions with stride [1 1] and padding [0 0 0 0] 65 'relu_conv10' ReLU ReLU 66 'pool10' Global Average Pooling Global average pooling 67 'prob' Softmax softmax 68 'ClassificationLayer_predictions' Classification Output crossentropyex with 'tench' and 999 other classes

You see that SqueezeNet consists of 68 layers. Like all DCNNs, SqueezeNet cascades convolutional operators followed by nonlinearities and pooling, or averaging. SqueezeNet expects an image input of size 227-by-227-by-3, which you can see with the following code.

snet.Layers(1)

ans = ImageInputLayer with properties: Name: 'data' InputSize: [227 227 3] Hyperparameters DataAugmentation: 'none' Normalization: 'zerocenter' NormalizationDimension: 'auto' Mean: [1×1×3 single]

Additionally, SqueezeNet is configured to recognized 1,000 different classes, which you can see with the following code.

snet.Layers(68)

ans = ClassificationOutputLayer with properties: Name: 'ClassificationLayer_predictions' Classes: [1000×1 categorical] OutputSize: 1000 Hyperparameters LossFunction: 'crossentropyex'

In a subsequent section, we will modify select layers of SqueezeNet in order to apply it to our classification problem.

SqueezeNet is designed to discriminate differences in images and classify the results. Therefore, in order to use SqueezeNet to classify radar returns, we must transform the 1-D radar return time series into an image. A common way to do this is to use a time-frequency representation (TFR). There are a number of choices for a time-frequency representation of a signal and which one is most appropriate depends on the signal characteristics. To determine which TFR may be appropriate for this problem, randomly choose and plot a few radar returns from each class.

```
rng default;
idxCylinder = randperm(50,2);
idxCone = randperm(50,2)+50;
```

It is evident that the radar returns previously shown are characterized by slowing varying changes punctuated by large transient decreases as described earlier. A wavelet transform is ideally suited to sparsely representing such signals. Wavelets shrink to localize transient phenomena with high temporal resolution and stretch to capture slowly varying signal structure. Obtain and plot the continuous wavelet transform of one of the cylinder returns.

`cwt(RCSReturns{:,idxCylinder(1)},'VoicesPerOctave',8)`

The CWT simultaneously captures both the slowly varying (low frequency) fluctuations and the transient phenomena. Contrast the CWT of the cylinder return with one from a cone target.

`cwt(RCSReturns{:,idxCone(2)},'VoicesPerOctave',8);`

Because of the apparent importance of the transients in determining whether the target return originates from a cylinder or cone target, we select the CWT as the ideal TFR to use. After obtaining the CWT for each target return, we make images from the CWT of each radar return. These images are resized to be compatible with SqueezeNet's input layer and we leverage SqueezeNet to classify the resulting images.

The helper function, `helpergenWaveletTFImg`

, obtains the CWT for each radar return, reshapes the CWT to be compatible with SqueezeNet, and writes the CWT as a jpeg file. To run `helpergenWaveletTFImg`

, choose a `parentDir`

where you have write permission. This example uses `tempdir`

, but you may use any folder on your machine where you have write permission. The helper function creates `Training`

and `Test`

set folders under `parentDir`

as well as creating `Cylinder`

and `Cone`

subfolders under both `Training`

and `Test`

. These folders are populated with jpeg images to be used as inputs to SqueezeNet.

parentDir = tempdir; helpergenWaveletTFImg(parentDir,RCSReturns,RCSReturnsTest)

Generating Time-Frequency Representations...Please Wait Creating Cylinder Time-Frequency Representations ... Done Creating Cone Time-Frequency Representations ... Done Creating Cylinder Time-Frequency Representations ... Done Creating Cone Time-Frequency Representations ... Done

Now use `imageDataStore`

to manage file access from the folders in order to train SqueezeNet. Create datastores for both the training and test data.

trainingData= imageDatastore(fullfile(parentDir,'Training'), 'IncludeSubfolders', true,... 'LabelSource', 'foldernames'); testData = imageDatastore(fullfile(parentDir,'Test'),'IncludeSubfolders',true,... 'LabelSource','foldernames');

In order to use SqueezeNet with this binary classification problem, we need to modify a couple layers. First, we change the last learnable layer in SqueezeNet (layer 64) to have the same number of 1-by-1 convolutions as our new number of classes, 2.

lgraphSqueeze = layerGraph(snet); convLayer = lgraphSqueeze.Layers(64); numClasses = numel(categories(trainingData.Labels)); newLearnableLayer = convolution2dLayer(1,numClasses, ... 'Name','binaryconv', ... 'WeightLearnRateFactor',10, ... 'BiasLearnRateFactor',10); lgraphSqueeze = replaceLayer(lgraphSqueeze,convLayer.Name,newLearnableLayer); classLayer = lgraphSqueeze.Layers(end); newClassLayer = classificationLayer('Name','binary'); lgraphSqueeze = replaceLayer(lgraphSqueeze,classLayer.Name,newClassLayer);

Finally, set the options for re-training SqueezeNet. Set the initial learn rate to 1e-4, set the maximum number of epochs to 15, and the minibatch size to 10. Use stochastic gradient descent with momentum.

ilr = 1e-4; mxEpochs = 15; mbSize =10; opts = trainingOptions('sgdm', 'InitialLearnRate', ilr, ... 'MaxEpochs',mxEpochs , 'MiniBatchSize',mbSize, ... 'Plots', 'training-progress','ExecutionEnvironment','cpu');

Train the network. If you have a compatible GPU, `trainNetwork`

automatically utilizes the GPU and training should complete in less than one minute. If you do not have a compatible GPU, `trainNetwork`

utilizes the CPU and training should take around five minutes. Training times do vary based on a number of factors. In this case, the training takes place on a cpu by setting the `ExecutionEnvironment`

parameter to `cpu`

.

CWTnet = trainNetwork(trainingData,lgraphSqueeze,opts);

Initializing input data normalization.

|========================================================================================| | Epoch | Iteration | Time Elapsed | Mini-batch | Mini-batch | Base Learning | | | | (hh:mm:ss) | Accuracy | Loss | Rate | |========================================================================================|

| 1 | 1 | 00:00:01 | 60.00% | 2.6639 | 1.0000e-04 |

| 5 | 50 | 00:00:42 | 100.00% | 0.0001 | 1.0000e-04 |

| 10 | 100 | 00:01:20 | 100.00% | 0.0002 | 1.0000e-04 |

| 15 | 150 | 00:01:58 | 100.00% | 2.2276e-05 | 1.0000e-04 |

|========================================================================================|

Use the trained network to predict target returns in the held-out test set.

predictedLabels = classify(CWTnet,testData,'ExecutionEnvironment','cpu'); accuracy = sum(predictedLabels == testData.Labels)/50*100

accuracy = 100

Plot the confusion chart along with the precision and recall. In this case, 100% of the test samples are classified correctly.

figure('Units','normalized','Position',[0.2 0.2 0.5 0.5]); ccDCNN = confusionchart(testData.Labels,predictedLabels); ccDCNN.Title = 'Confusion Chart'; ccDCNN.ColumnSummary = 'column-normalized'; ccDCNN.RowSummary = 'row-normalized';

In the final section of this example, an LSTM workflow is described. First the LSTM layers are defined:

LSTMlayers = [ ... sequenceInputLayer(1) bilstmLayer(100,'OutputMode','last') fullyConnectedLayer(2) softmaxLayer classificationLayer ]; options = trainingOptions('adam', ... 'MaxEpochs',30, ... 'MiniBatchSize', 150, ... 'InitialLearnRate', 0.01, ... 'GradientThreshold', 1, ... 'plots','training-progress', ... 'Verbose',false,'ExecutionEnvironment','cpu'); trainLabels = repelem(categorical({'cylinder','cone'}),[50 50]); trainLabels = trainLabels(:); trainData = num2cell(table2array(RCSReturns)',2); testData = num2cell(table2array(RCSReturnsTest)',2); testLabels = repelem(categorical({'cylinder','cone'}),[25 25]); testLabels = testLabels(:); RNNnet = trainNetwork(trainData,trainLabels,LSTMlayers,options);

The accuracy for this system is also plotted.

predictedLabels = classify(RNNnet,testData,'ExecutionEnvironment','cpu'); accuracy = sum(predictedLabels == testLabels)/50*100

accuracy = 100

This example presents a workflow for performing radar target classification using machine and deep learning techniques. Although this example used synthesized data to do training and testing, it can be easily extended to accommodate real radar returns. Because of the signal characteristics, wavelet techniques were used for both the machine learning and CNN approaches.

With this dataset we were also obtained to achieve similar accuracy by just feeding the raw data into a LSTM. In more complicated datasets, the raw data may be too inherently variable for the model to learn robust features from the raw data and you may have to resort to feature extraction prior to using a LSTM.