www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/PhysiologicSignalAnalysisExample.m
%% Wavelet Analysis of Physiologic Signals % This example shows how to use wavelets to analyze physiologic signals. % % Physiologic signals are frequently nonstationary meaning that their % frequency content changes over time. In many applications, these changes % are the events of interest. % % Wavelets decompose signals into time-varying frequency (scale) % components. Because signal features are often localized in time and % frequency, analysis and estimation are easier when working with % sparser (reduced) representations. % % This example presents a few illustrative cases where the ability of % wavelets to provide a local time-frequency analysis of signals is % beneficial. %% R Wave Detection in the Electrocardiogram with MODWT % The QRS complex consists of three deflections in the electrocardiogram % (ECG) waveform. The QRS complex reflects the depolarization of the right % and left ventricles and is the most prominent feature of the human ECG. % % Load and plot an ECG waveform where the R peaks of the QRS complex have % been annotated by two or more cardiologists. The ECG data and annotations % are taken from the MIT-BIH Arrythmia Database. The data are sampled at % 360 Hz. load mit200; figure plot(tm,ecgsig) hold on plot(tm(ann),ecgsig(ann),'ro') xlabel('Seconds'); ylabel('Amplitude') title('Subject - MIT-BIH 200') %% % You can use wavelets to build an automatic QRS detector for use in % applications like R-R interval estimation. % % There are two keys for using wavelets as general feature detectors: % % * The wavelet transform separates signal components into different % frequency bands enabling a sparser representation of the signal. % % * You can often find a wavelet which resembles the feature you are trying % to detect. % % The 'sym4' wavelet resembles the QRS complex, which makes it a good % choice for QRS detection. To illustrate this more clearly, extract a QRS % complex and plot the result with a dilated and translated 'sym4' % wavelet for comparison. qrsEx = ecgsig(4560:4810); [mpdict,~,~,longs] = wmpdictionary(numel(qrsEx),'lstcpt',{{'sym4',3}}); figure plot(qrsEx) hold on plot(2*circshift(mpdict(:,11),[-2 0]),'r') axis tight legend('QRS Complex','Sym4 Wavelet') title('Comparison of Sym4 Wavelet and QRS Complex') %% % Use the maximal overlap discrete wavelet transform (MODWT) to enhance the % R peaks in the ECG waveform. The MODWT is an undecimated wavelet % transform, which handles arbitrary sample sizes. % % First, decompose the ECG waveform down to level 5 using the default % 'sym4' wavelet. Then, reconstruct a frequency-localized version of the % ECG waveform using only the wavelet coefficients at scales 4 and 5. The % scales correspond to the following approximate frequency bands. % % * Scale 4 -- [11.25, 22.5) Hz % * Scale 5 -- [5.625, 11.25) Hz. % % This covers the passband shown to maximize QRS energy. wt = modwt(ecgsig,5); wtrec = zeros(size(wt)); wtrec(4:5,:) = wt(4:5,:); y = imodwt(wtrec,'sym4'); %% % Use the squared absolute values of the signal approximation built from % the wavelet coefficients and employ a peak finding algorithm to identify % the R peaks. % % If you have the Signal Processsing Toolbox(TM), you can use % |findpeaks| to locate the peaks. Plot the R-peak waveform obtained % with the wavelet transform annotated with the automatically-detected peak % locations. y = abs(y).^2; [qrspeaks,locs] = findpeaks(y,tm,'MinPeakHeight',0.35,... 'MinPeakDistance',0.150); figure; plot(tm,y) hold on plot(locs,qrspeaks,'ro') xlabel('Seconds') title('R Peaks Localized by Wavelet Transform with Automatic Annotations') %% % Add the expert annotations to the R-peak waveform. Automatic peak % detection times are considered accurate if within 150 msec of the true % peak ($\pm 75$ msec). plot(tm(ann),y(ann),'k*') title('R peaks Localized by Wavelet Transform with Expert Annotations') %% % At the command line, you can compare the values of |tm(ann)| and |locs|, % which are the expert times and automatic peak detection times % respectively. Enhancing the R peaks with the wavelet transform results in % a hit rate of 100% and no false positives. The calculated heart rate % using the wavelet transform is 88.60 beats/minute compared to 88.72 % beats/minute for the annotated waveform. %% % If you try to work on the square magnitudes of the original data, you % find the capability of the wavelet transform to isolate the R peaks makes % the detection problem much easier. Working on the raw data can cause % misidentifications such as when the squared S-wave peak exceeds the % R-wave peak around 10.4 seconds. figure plot(tm,ecgsig,'k--') hold on plot(tm,y,'r','linewidth',1.5) plot(tm,abs(ecgsig).^2,'b') plot(tm(ann),ecgsig(ann),'ro','markerfacecolor',[1 0 0]) set(gca,'xlim',[10.2 12]) legend('Raw Data','Wavelet Reconstruction','Raw Data Squared', ... 'Location','SouthEast'); xlabel('Seconds') %% % Using |findpeaks| on the squared magnitudes of the raw data results in % twelve false positives. [qrspeaks,locs] = findpeaks(ecgsig.^2,tm,'MinPeakHeight',0.35,... 'MinPeakDistance',0.150); %% % In addition to switches in polarity of the R peaks, the ECG is often % corrupted by noise. load mit203; figure plot(tm,ecgsig) hold on plot(tm(ann),ecgsig(ann),'ro') xlabel('Seconds'); ylabel('Amplitude') title('Subject - MIT-BIH 203 with Expert Annotations') %% % Use the MODWT to isolate the R peaks. Use |findpeaks| to determine the % peak locations. Plot the R-peak waveform along with the expert and % automatic annotations. wt = modwt(ecgsig,5); wtrec = zeros(size(wt)); wtrec(4:5,:) = wt(4:5,:); y = imodwt(wtrec,'sym4'); y = abs(y).^2; [qrspeaks,locs] = findpeaks(y,tm,'MinPeakHeight',0.1,... 'MinPeakDistance',0.150); figure plot(tm,y) title('R-Waves Localized by Wavelet Transform') hold on hwav = plot(locs,qrspeaks,'ro'); hexp = plot(tm(ann),y(ann),'k*'); xlabel('Seconds') legend([hwav hexp],'Automatic','Expert','Location','NorthEast'); %% % The hit rate is again 100% with zero false alarms. %% % The previous examples used a very simple wavelet QRS detector based on a % signal approximation constructed from <matlab:doc('modwt') modwt>. The % goal was to demonstrate the ability of the wavelet transform to isolate % signal components, not to build the most robust wavelet-transform-based % QRS detector. It is possible, for example, to exploit the fact that the % wavelet transform provides a multiscale analysis of the signal to enhance % peak detection. Examine the scale 4 and 5 magnitude-squared wavelet % details plotted along with R peak times as annotated by the experts. The % level-4 details are shifted for visualization. ecgmra = modwtmra(wt); figure plot(tm,ecgmra(5,:).^2,'b'); hold on; plot(tm,ecgmra(4,:).^2+0.6,'b'); set(gca,'xlim',[14.3 25.5]); timemarks = repelem(tm(ann),2); N = numel(timemarks); markerlines = reshape(repmat([0;1],1,N/2),N,1); h = stem(timemarks,markerlines,'k--'); h.Marker = 'none'; set(gca,'ytick',[0.1 0.6]); set(gca,'yticklabels',{'D5','D4'}) xlabel('Seconds') title('Magnitude-Squared Level 4 and 5 Details') %% % You see that peaks in the level 4 and level 5 details tend to co-occur. A % more advanced wavelet peak-finding algorithm could exploit this by using % information from multiple scales simultaneously. %% Time-Varying Wavelet Coherence Analysis of Brain Dynamics % Fourier-domain coherence is a well-established technique for measuring % the linear correlation between two stationary processes as a function of % frequency on a scale from 0 to 1. Because wavelets provide local % information about data in time and scale (frequency), wavelet-based % coherence allows you to measure time-varying correlation as a function of % frequency. In other words, a coherence measure suitable for nonstationary % processes. % % To illustrate this, examine near-infrared spectroscopy (NIRS) data % obtained in two human subjects. NIRS measures brain activity by % exploiting the different absorption characteristics of oxygenated and % deoxygenated hemoglobin. The data is taken from Cui, Bryant, & Reiss % (2012) and was kindly provided by the authors for this example. The % recording site was the superior frontal cortex for both subjects. The % data is sampled at 10 Hz. % % In the experiment, the subjects alternatively cooperated and % competed on a task. The period of the task was seven seconds. load NIRSData; figure plot(tm,[NIRSData(:,1) NIRSData(:,2)]) legend('Subject 1','Subject 2','Location','NorthWest') xlabel('Seconds') title('NIRS Data') grid on; %% % Examining the time-domain data, it is not clear what oscillations are % present in the individual time series, or what oscillations are common to % both data sets. Use wavelet analysis to answer both questions. cwt(NIRSData(:,1),1/10,'bump'); figure; cwt(NIRSData(:,2),1/10,'bump'); %% % The CWT analyses reveal strong frequency-modulated oscillations in both % datasets around 1 Hz. These are due to the cardiac cycles in the two % subjects. Additionally, there appears to be a weaker oscillation in both % datasets around 0.15 Hz. This activity is stronger and more consistent in % subject 1 than subject 2. Wavelet coherence can enhance the detection of % weak oscillations that are jointly present in two time series. [wcoh,~,F] = wcoherence(NIRSData(:,1),NIRSData(:,2),10); figure; surf(tm,F,abs(wcoh).^2); view(0,90); shading interp; axis tight; hc = colorbar; hc.Label.String = 'Coherence'; title('Wavelet Coherence') xlabel('Seconds'),ylabel('Hz'); ylim([0 2.5]); set(gca,'ytick',[0.15 1.2 2]) %% % In the wavelet coherence, there is a strong correlation around % 0.15 Hz. This is in the frequency band corresponding to the experimental % task and represents task-related coherent oscillations in brain activity % in the two subjects. Add to the plot time markers which indicate % two task periods. The period between the tasks is a rest period. taskbd = [245 1702 2065 3474]; tvec = repelem(tm(taskbd),2); yvec = [0 max(F)]'; yvec = reshape(repmat(yvec,1,4),8,1); hold on; stemPlot = stem(tvec,yvec,'w--','linewidth',2); stemPlot.Marker = 'none'; %% % This example used <matlab:doc('cwt') cwt> to obtain and plot the a % time-frequency analysis of the individual NIRS time series. The example % also used <matlab:doc('wcoherence') wcoherence> to obtain the wavelet % coherence of the two time series. The use of wavelet coherence often % enables you to detect coherent oscillatory behavior in two time series % which may be fairly weak in each individual series. Consult Cui, Bryant, % & Reiss (2012) for a more detailed wavelet-coherence analysis of this % data. %% Time-Frequency Analysis of Otoacoustic Emission Data with the CWT % Otoacoustic emissions (OAEs) are narrowband oscillatory signals emitted % by the cochlea (inner ear) and their presence is indicative of normal % hearing. Load and plot some example OAE data. The data are sampled at 20 % kHz. load dpoae figure plot(t.*1000,dpoaets) xlabel('Milliseconds') ylabel('Amplitude') %% % The emission was evoked by a stimulus beginning at 25 milliseconds and % ending at 175 milliseconds. Based on the experimental parameters, the % emission frequency should be 1230 Hz. Obtain and plot the CWT as a % function of time and frequency. Use the default analytic Morse wavelet % with 16 voices per octave. [dpoaeCWT,f] = cwt(dpoaets,2e4,'VoicesPerOctave',16); helperCWTTimeFreqPlot(dpoaeCWT,t.*1000,f,... 'surf','CWT of OAE','milliseconds','Hz') %% % You can investigate the time evolution of the OAE by finding the CWT % coefficients closest in frequency to 1230 Hz and examining their % magnitudes as a function of time. Plot the magnitudes along with time % markers designating the beginning and end of the evoking stimulus. [~,idx1230] = min(abs(f-1230)); cfsOAE = dpoaeCWT(idx1230,:); plot(t.*1000,abs(cfsOAE)); hold on AX = gca; plot([25 25],[AX.YLim(1) AX.YLim(2)],'r') plot([175 175],[AX.YLim(1) AX.YLim(2)],'r') xlabel('msec') title('CWT Coefficient Magnitudes') %% % There is some delay between the onset of the evoking stimulus and the % OAE. Once the evoking stimulus is terminated, the OAE immediately begins % to decay in magnitude. %% % Another way to isolate the emission is to use the inverse CWT to % reconstruct a frequency-localized approximation in the time domain. % % Construct a frequency-localized emission approximation by extracting the % CWT coefficients corresponding to frequencies between 1150 and 1350 Hz. % Use these coefficients and invert the CWT. Plot the original data along % with the reconstruction and markers indicating the beginning and end of % the evoking stimulus. frange = [1150 1350]; xrec = icwt(dpoaeCWT,f,frange); figure plot(t.*1000,dpoaets) hold on; plot(t.*1000,xrec,'r') AX = gca; ylimits = AX.YLim; plot([25 25],ylimits,'k') plot([175 175],ylimits,'k') grid on; xlabel('Milliseconds') ylabel('Amplitude') title('Frequency-Localized Reconstruction of Emission') %% % In the time-domain data, you clearly see how the emission ramps on and % off at the application and termination of the evoking stimulus. % % It is important to note that even though a range of frequencies were % selected for the reconstruction, the analytic wavelet transform actually % encodes the exact frequency of the emission. To demonstrate this, take % the Fourier transform of the emission approximation reconstructed from % the analytic CWT. xdft = fft(xrec); freq = 0:2e4/numel(xrec):1e4; xdft = xdft(1:numel(xrec)/2+1); figure plot(freq,abs(xdft)); xlabel('Hz'); ylabel('Magnitude') title('Fourier Transform of CWT-Based Signal Approximation'); [~,maxidx] = max(abs(xdft)); fprintf('The frequency is %4.2f Hz\n',freq(maxidx)); %% % This example used |cwt| to obtain a time-frequency analysis of the OAE % data and <matlab:doc('icwt') icwt> to obtain a frequency-localized % approximation to the signal. %% References % % Cui, X., Bryant, D.M., and Reiss. A.L. "NIRS-Based hyperscanning reveals % increased interpersonal coherence in superior frontal cortex during % cooperation", Neuroimage, 59(3), 2430-2437, 2012. % % Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, % Mietus JE, Moody GB, Peng C-K, Stanley HE. "PhysioBank, PhysioToolkit, % and PhysioNet: Components of a New Research Resource for Complex % Physiologic Signals." Circulation 101(23):e215-e220, 2000. % |http://circ.ahajournals.org/cgi/content/full/101/23/e215| % % Mallat, S. "A Wavelet Tour of Signal Processing: The Sparse Way", % Academic Press, 2009. % % Moody, G.B. "Evaluating ECG Analyzers". % |http://www.physionet.org/physiotools/wfdb/doc/wag-src/eval0.tex| % % Moody GB, Mark RG. "The impact of the MIT-BIH Arrhythmia Database." IEEE % Eng in Med and Biol 20(3):45-50 (May-June 2001). %% Appendix % The following helper functions are used in this example. % % * <matlab:edit('helperCWTTimeFreqPlot.m') helperCWTTimeFreqPlot.m> % * <matlab:edit('helperCWTTimeFreqVector.m') helperCWTTimeFreqVector.m>