Data Description

This example uses ECG data obtained from three groups, or classes, of people: persons with cardiac arrhythmia, persons with congestive heart failure, and persons with normal sinus rhythms. The example uses 162 ECG recordings from three PhysioNet databases:MIT-BIH Arrhythmia Database[3][5],MIT-BIH Normal Sinus Rhythm Database[3], andThe BIDMC Congestive Heart Failure Database[2][3]. In total, there are 96 recordings from persons with arrhythmia, 30 recordings from persons with congestive heart failure, and 36 recordings from persons with normal sinus rhythms. The goal is to train a classifier to distinguish between arrhythmia (ARR), congestive heart failure (CHF), and normal sinus rhythm (NSR).

Download Data

The first step is to download the data from theGitHub回购itory。To download the data, clickCode和selectDownload ZIP。Save the filephysionet_ECG_data-master.zipin a folder where you have write permission. The instructions for this example assume you have downloaded the file to your temporary directory, (tempdirin MATLAB). Modify the subsequent instructions for unzipping and loading the data if you choose to download the data in folder different fromtempdir。

The filephysionet_ECG_data-master.zipcontains

和ECGData.zip contains

ECGData.mat holds the data used in this example. The .txt file, Modified_physionet_data.txt, is required by PhysioNet's copying policy and provides the source attributions for the data as well as a description of the pre-processing steps applied to each ECG recording.

Load Files

If you followed the download instructions in the previous section, enter the following commands to unzip the two archive files.

unzip(fullfile(tempdir,'physionet_ECG_data-master.zip'),tempdir) unzip(fullfile(tempdir,'physionet_ECG_data-master','ECGData.zip'),。..fullfile(tempdir,'ECGData'))

After you unzip the ECGData.zip file, load the data into MATLAB.

load(fullfile(tempdir,'ECGData','ECGData.mat'))

ECGDatais a structure array with two fields:Data和Labels。Datais a 162-by-65536 matrix where each row is an ECG recording sampled at 128 hertz. Each ECG time series has a total duration of 512 seconds.Labelsis a 162-by-1 cell array of diagnostic labels, one for each row ofData。The three diagnostic categories are: 'ARR' (arrhythmia), 'CHF' (congestive heart failure), and 'NSR' (normal sinus rhythm).

Create Training and Test Data

Randomly split the data into two sets - training and test data sets. The helper functionhelperRandomSplitperforms the random split.helperRandomSplitaccepts the desired split percentage for the training data andECGData。ThehelperRandomSplitfunction outputs two data sets along with a set of labels for each. Each row oftrainData和testDatais an ECG signal. Each element oftrainLabels和testLabelscontains the class label for the corresponding row of the data matrices. In this example, we randomly assign 70% percent of the data in each class to the training set. The remaining 30% is held out for testing (prediction) and are assigned to the test set.

percent_train = 70; [trainData,testData,trainLabels,testLabels] =。..helperRandomSplit(percent_train,ECGData);

There are 113 records in thetrainDataset and 49 records intestData。通过设计训练数据包含69.75% (113/162) of the data. Recall that the ARR class represents 59.26% of the data (96/162), the CHF class represents 18.52% (30/162), and the NSR class represents 22.22% (36/162). Examine the percentage of each class in the training and test sets. The percentages in each are consistent with the overall class percentages in the data set.

Ctrain = countcats(categorical(trainLabels))./numel(trainLabels).*100 Ctest = countcats(categorical(testLabels))./numel(testLabels).*100

Ctrain = 59.2920 18.5841 22.1239 Ctest = 59.1837 18.3673 22.4490

Plot Samples

Plot the first few thousand samples of four randomly selected records fromECGData。The helper functionhelperPlotRandomRecordsdoes this.helperPlotRandomRecordsacceptsECGData和a random seed as input. The initial seed is set at 14 so that at least one record from each class is plotted. You can executehelperPlotRandomRecordswithECGDataas the only input argument as many times as you wish to get a sense of the variety of ECG waveforms associated with each class. You can find the source code for this and all helper functions in the Supporting Functions section at the end of this example.

helperPlotRandomRecords(ECGData,14)

Wavelet Time Scattering

The key parameters to specify in a wavelet time scattering network 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 network 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 150 seconds.

N = size(ECGData.Data,2); sn = waveletScattering('SignalLength',N,'InvarianceScale',150,'SamplingFrequency',128);

You can visualize the wavelet filters in the two filter banks with the following.

[fb,f,filterparams] = filterbank(sn); subplot(211) plot(f,fb{2}.psift) xlim([0 128]) gridontitle('1st Filter Bank Wavelet Filters'); subplot(212) plot(f,fb{3}.psift) xlim([0 128]) gridontitle('2nd Filter Bank Wavelet Filters'); xlabel('Hz');

To demonstrate the invariance scale, obtain the inverse Fourier transform of the scaling function and center it at 0 seconds in time. The two vertical black lines mark the -75 and 75 second boundaries. Additionally, plot the real and imaginary parts of the coarsest-scale (lowest frequency) wavelet from the first filter bank. Note that the coarsest-scale wavelet does not exceed the invariant scale determined by the time support of the scaling function. This is an important property of wavelet time scattering.

figure; phi = ifftshift(ifft(fb{1}.phift)); psiL1 = ifftshift(ifft(fb{2}.psift(:,end))); t = (-2^15:2^15-1).*1/128; scalplt = plot(t,phi); holdongridonylim([-1.5e-4 1.6e-4]); plot([-75 -75],[-1.5e-4 1.6e-4],'k--'); plot([75 75],[-1.5e-4 1.6e-4],'k--'); xlabel('Seconds'); ylabel('Amplitude'); wavplt = plot(t,[real(psiL1) imag(psiL1)]); legend([scalplt wavplt(1) wavplt(2)],{'Scaling Function','Wavelet-Real Part','Wavelet-Imaginary Part'}); title({'Scaling Function';'Coarsest-Scale Wavelet First Filter Bank'}) holdoff

After constructing the scattering network, obtain the scattering coefficients for the training data as a matrix. When you runfeatureMatrixwith multiple signals, each column is treated a single signal.

scat_features_train = featureMatrix(sn,trainData');

The output offeatureMatrixin this case is 409-by-16-by-113. Each page of the tensor,scat_features_train, is the scattering transform of one signal. The wavelet scattering transform is critically downsampled in time based on the bandwidth of the scaling function. In this case, this results in 16 time windows for each of the 409 scattering paths.

In order to obtain a matrix compatible with the SVM classifier, reshape the multisignal scattering transform into a matrix where each column corresponds to a scattering path and each row is a scattering time window. In this case, you obtain 1808 rows because there are 16 time windows for each of the 113 signals in the training data.

Nwin = size(scat_features_train,2); scat_features_train = permute(scat_features_train,[2 3 1]); scat_features_train = reshape(scat_features_train,。..size(scat_features_train,1)*size(scat_features_train,2),[]);

Repeat the process for the test data. Initially,scat_features_testis 409-by-16-by-49 because there are 49 ECG waveforms in the training set. After reshaping for the SVM classifier, the feature matrix is 784-by-416.

scat_features_test = featureMatrix(sn,testData'); scat_features_test = permute(scat_features_test,[2 3 1]); scat_features_test = reshape(scat_features_test,。..size(scat_features_test,1)*size(scat_features_test,2),[]);

Because for each signal we obtained 16 scattering windows, we need to create labels to match the number of windows. The helper functioncreateSequenceLabelsdoes this based on the number of windows.

[sequence_labels_train,sequence_labels_test] = createSequenceLabels(Nwin,trainLabels,testLabels);

Cross Validation

For classification, two analyses are performed. The first uses all the scattering data and fits a multi-class SVM with a quadratic kernel. In total there 2592 scattering sequences in the entire dataset, 16 for each of the 162 signals. The error rate, or loss, is estimated using 5-fold cross validation.

scat_features = [scat_features_train; scat_features_test]; allLabels_scat = [sequence_labels_train; sequence_labels_test]; rng(1); template = templateSVM(。..'KernelFunction','polynomial',。..'PolynomialOrder', 2,。..'KernelScale','auto',。..'BoxConstraint', 1,。..'Standardize', true); classificationSVM = fitcecoc(。..scat_features,。..allLabels_scat,。..“学习者”, template,。..'Coding','onevsone',。..'ClassNames', {'ARR';'CHF';'NSR'}); kfoldmodel = crossval(classificationSVM,'KFold', 5);

Compute the loss and the confusion matrix. Display the accuracy.

predLabels = kfoldPredict(kfoldmodel); loss = kfoldLoss(kfoldmodel)*100; confmatCV = confusionmat(allLabels_scat,predLabels) fprintf('Accuracy is %2.2f percent.\n',100-loss);

confmatCV = 1535 0 1 2 478 0 0 0 576 Accuracy is 99.88 percent.

The accuracy is 99.88%, which is quite good but the actual results are better considering that here each time window is classified separately. There are 16 separate classifications for each signal. Use a simple majority vote to obtain a single class prediction for each scattering representation.

classes = categorical({'ARR','CHF','NSR'}); [ClassVotes,ClassCounts] = helperMajorityVote(predLabels,[trainLabels; testLabels],classes);

Determine the actual cross-validation accuracy based on the mode of the class predictions for each set of scattering time windows. If the mode is not unique for a given set,helperMajorityVotereturns a classification error indicated by'NoUniqueMode'。This results in an extra column in the confusion matrix, which in this case is all zeros because a unique mode exists for each set of scattering predictions.

CVaccuracy = sum(eq(ClassVotes,categorical([trainLabels; testLabels])))/162*100; fprintf('True cross-validation accuracy is %2.2f percent.\n',CVaccuracy); MVconfmatCV = confusionmat(categorical([trainLabels; testLabels]),ClassVotes); MVconfmatCV

True cross-validation accuracy is 100.00 percent. MVconfmatCV = 96 0 0 0 0 30 0 0 0 0 36 0 0 0 0 0

Scattering has correctly classified all the signals in the cross-validated model. If you examineClassCounts, you see that the 2 misclassified time windows inconfmatCVare attributable to 2 signals where 15 of the 16 scattering windows were classified correctly.

SVM Classification

For the next analysis, we fit a multi-class quadratic SVM to the training data only (70%) and then use that model to make predictions on the 30% of the data held out for testing. There are 49 data records in the test set. Use a majority vote on the individual scattering windows.

model = fitcecoc(。..scat_features_train,。..sequence_labels_train,。..“学习者”, template,。..'Coding','onevsone',。..'ClassNames', {'ARR','CHF','NSR'}); predLabels = predict(model,scat_features_test); [TestVotes,TestCounts] = helperMajorityVote(predLabels,testLabels,classes); testaccuracy = sum(eq(TestVotes,categorical(testLabels)))/numel(testLabels)*100; fprintf(的测试精度是2.2 f % %。\ n',testaccuracy); confusionchart(categorical(testLabels),TestVotes)

The test accuracy is 97.96 percent.

The classification accuracy on the test dataset is approximately 98%. The confusion matrix shows that one CHF record is misclassified as ARR. All 48 other signals are correctly classified.

Summary

This example used wavelet time scattering and an SVM classifier to classify ECG waveforms into one of three diagnostic classes. Wavelet scattering proved to be a powerful feature extractor, which required only a minimal set of user-specified parameters to yield a set of robust features for classification. Compare this with the exampleSignal Classification Using Wavelet-Based Features and Support Vector Machines, which required a significant amount of expertise to handcraft features to use in classification. With wavelet time scattering, you are only required to specify the scale of the time invariance, the number of filter banks (or wavelet transforms), and the number of wavelets per octave. The combination of a wavelet scattering transform and an SVM classifier yielded 100% classification on a cross-validated model and 98% correct classification when applying an SVM to the scattering transforms of a hold-out test set.

References

Supporting Functions

helperPlotRandomRecordsPlots four ECG signals randomly chosen fromECGData。

functionhelperPlotRandomRecords(ECGData,randomSeed)% This function is only intended to support the XpwWaveletMLExample. It may% change or be removed in a future release.ifnargin==2 rng(randomSeed)endM = size(ECGData.Data,1); idxsel = randperm(M,4);fornumplot = 1:4 subplot(2,2,numplot) plot(ECGData.Data(idxsel(numplot),1:3000)) ylabel('Volts')ifnumplot > 2 xlabel('Samples')endtitle(ECGData.Labels{idxsel(numplot)})endend

helperMajorityVoteFinds the mode in the predicted class labels for each set of scattering time windows. The function returns both the class label modes and the number of class predictions for each set of scattering time windows. If there is no unique mode,helperMajorityVotereturns a class label of "error" to indicate that set of scattering windows is a classification error.

function[ClassVotes,ClassCounts] = helperMajorityVote(predLabels,origLabels,classes)% This function is in support of ECGWaveletTimeScatteringExample. It may% change or be removed in a future release.% Make categorical arrays if the labels are not already categoricalpredLabels = categorical(predLabels); origLabels = categorical(origLabels);% Expects both predLabels and origLabels to be categorical vectorsNpred = numel(predLabels); Norig = numel(origLabels); Nwin = Npred/Norig; predLabels = reshape(predLabels,Nwin,Norig); ClassCounts = countcats(predLabels); [mxcount,idx] = max(ClassCounts); ClassVotes = classes(idx);% Check for any ties in the maximum values and ensure they are marked as% error if the mode occurs more than oncemodecnt = modecount(ClassCounts,mxcount); ClassVotes(modecnt>1) = categorical({'error'}); ClassVotes = ClassVotes(:);%-------------------------------------------------------------------------functionmodecnt = modecount(ClassCounts,mxcount) modecnt = Inf(size(ClassCounts,2),1);fornc = 1:size(ClassCounts,2) modecnt(nc) = histc(ClassCounts(:,nc),mxcount(nc));endend% EOFend