www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/wcoherdemo.m
%% Compare Time-Frequency Content in Signals with Wavelet Coherence % This example shows how to use wavelet coherence and the wavelet % cross-spectrum to identify time-localized common oscillatory behavior in % two time series. The example also compares the wavelet coherence and % cross-spectrum against their Fourier counterparts. You must have the % Signal Processing Toolbox(TM) to run the examples using |mscohere| and % |cpsd|. % % Many applications involve identifying and characterizing common patterns % in two time series. In some situations, common behavior in two time % series results from one time series driving or influencing the other. In % other situations, the common patterns result from some unobserved % mechanism influencing both time series. % % For jointly stationary time series, the standard techniques for % characterizing correlated behavior in time or frequency are % cross-correlation, the (Fourier) cross-spectrum, and coherence. However, % many time series are nonstationary, meaning that their frequency content % changes over time. For these time series, it is important to have a % measure of correlation or coherence in the time-frequency plane. % % You can use wavelet coherence to detect common time-localized % oscillations in nonstationary signals. In situations where it is natural % to view one time series as influencing another, you can use the phase of % the wavelet cross-spectrum to identify the relative lag between the two % time series. %% Locate Common Time-Localized Oscillations and Determine Phase Lag % For the first example, use two signals consisting of time-localized % oscillations at 10 and 75 Hz. The signals are six seconds in duration and % are sampled at 1000 Hz. The 10-Hz oscillation in the two signals overlaps % between 1.2 and 3 seconds. The overlap for the 75-Hz oscillation occurs % between 0.4 and 4.4 seconds. The 10 and 75-Hz components are delayed 1/4 % of a cycle in the Y-signal. This means there is a $\pi/2$ (90 degree) % phase lag between the 10 and 75-Hz components in the two signals. Both % signals are corrupted by additive white Gaussian noise. %% load wcoherdemosig1; subplot(2,1,1) plot(t,x1); title('X Signal') grid on; ylabel('Amplitude'); subplot(2,1,2) plot(t,y1); title('Y Signal'); ylabel('Amplitude'); grid on; xlabel('Seconds'); %% % Obtain the wavelet coherence and plot the result. Enter the sampling % frequency (1000 Hz) to obtain a time-frequency plot of the wavelet % coherence. In regions of the time-frequency plane where coherence exceeds % 0.5, the phase from the wavelet cross-spectrum is used to indicate the % relative lag between coherent components. Phase is indicated by arrows % oriented in a particular direction. Note that a 1/4 cycle lag in the % Y-signal at a particular frequency is indicated by an arrow pointing % vertically. figure; wcoherence(x1,y1,1000); %% % The two time-localized regions of coherent oscillatory behavior at 10 and % 75 Hz are evident in the plot of the wavelet coherence. The phase % relationship is shown by the orientation of the arrows in the regions of % high coherence. In this example, you see that the wavelet cross-spectrum % captures the $\pi/2$ (1/4 cycle) phase lag between the two signals at 10 % and 75 Hz. % % The white dashed line shows the cone of influence where edge effects % become significant at different frequencies (scales). Areas of high % coherence occuring outside or overlapping the cone of influence should be % interpreted with caution. %% % Repeat the same analysis using the Fourier magnitude-squared coherence % and cross-spectrum. figure; mscohere(x1,y1,500,450,[],1000); [Pxy,F] = cpsd(x1,y1,500,450,[],1000); Phase = (180/pi)*angle(Pxy); figure; plot(F,Phase); title('Cross-Spectrum Phase'); xlabel('Hz'); ylabel('Phase (Degrees)'); grid on; ylim([0 180]); xlim([0 200]); hold on; plot([10 10],[0 180],'r--'); plot([75 75],[0 180],'r--'); plot(F,90*ones(size(F)),'r--'); hold off; %% % The Fourier magnitude-squared coherence obtained from |mscohere| clearly % identifies the coherent oscillations at 10 and 75 Hz. In the phase plot % of the Fourier cross-spectrum, the vertical red dashed lines mark 10 and % 75 Hz while the horizontal line marks an angle of 90 degrees. You see % that the phase of the cross-spectrum does a reasonable job of capturing % the relative phase lag between the components. % % However, the time-dependent nature of the coherent behavior is completely % obscured by these techniques. For nonstationary signals, characterizing % coherent behavior in the time-frequency plane is much more informative. %% % The following example repeats the preceding one while changing the phase % relationship between the two signals. In this case, the 10-Hz component % in the Y-signal is delayed by 3/8 of a cycle ($\frac{3 \pi}{4}$ radians). % The 75-Hz component in the Y-signal is delayed by 1/8 of a cycle % ($\frac{\pi}{4}$). Plot the wavelet coherence and threshold the phase % display to only show areas where the coherence exceeds 0.75. load wcoherdemosig2; wcoherence(x2,y2,1000,'phasedisplaythreshold',0.75); %% % Phase estimates obtained from the wavelet cross-spectrum capture the % the relative lag between the two time series at 10 and 75 Hz. %% % If you prefer to view data in terms of periods rather than frequency, you % can use a MATLAB duration object to provide |wcoherence| with a sample % time. wcoherence(x2,y2,seconds(0.001),'phasedisplaythreshold',0.75); %% % Note that cone of influence has inverted because the wavelet coherence is % now plotted in terms of periods. %% Determine Coherent Oscillations and Delay in Climate Data % Load and plot the El Nino Region 3 data and deasonalized All Indian % Rainfall Index from 1871 to late 2003. The data are sampled monthly. The % Nino 3 time series is a record of monthly sea surface temperature % anomalies in degrees Celsius recorded from the equatorial Pacific at % longitudes 90 degrees west to 150 degrees west and latitudes 5 degrees % north to 5 degrees south. The All Indian Rainfall Index represents % average Indian rainfalls in millimeters with seasonal components removed. load ninoairdata; figure; subplot(2,1,1) plot(datayear,nino); title('El Nino Region 3 -- SST Anomalies'); ylabel('Degrees Celsius'); axis tight; subplot(2,1,2); plot(datayear,air); axis tight; title('Deasonalized All Indian Rainfall Index'); ylabel('Millimeters'); xlabel('Year'); %% % Plot the wavelet coherence with phase estimates. For this data, it is % more natural to look at oscillations in terms of their periods in years. % Enter the sampling interval (period) as duration object with units of % years so that the output periods are in years. Show the relative lag % between the two climate time series only where the magnitude-squared % coherence exceeds 0.7. figure; wcoherence(nino,air,years(1/12),'phasedisplaythreshold',0.7); %% % The plot shows time-localized areas of strong coherence occuring in % periods that correspond to the typical El Nino cycles of 2 to 7 years. % The plot also shows that there is an approximate 3/8-to-1/2 cycle delay % between the two time series at those periods. This indicates that periods % of sea warming consistent with El Nino recorded off the coast of South % America are correlated with rainfall amounts in India approximately % 17,000 km away, but that this effect is delayed by approximately 1/2 a % cycle (1 to 3.5 years). %% Find Coherent Oscillations in Brain Activity % In the previous examples, it was natural to view one time series as % influencing the other. In these cases, examining the lead-lag % relationship between the data is informative. In other cases, it is more % natural to examine the coherence alone. % % For an example, consider 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 approximately 7.5 seconds. load NIRSData; figure plot(tm,NIRSData(:,1)) hold on plot(tm,NIRSData(:,2),'r') legend('Subject 1','Subject 2','Location','NorthWest') xlabel('Seconds') title('NIRS Data') grid on; hold off; %% % Obtain the wavelet coherence as a function of time and frequency. You can % use |wcoherence| to output the wavelet coherence, cross-spectrum, % scale-to-frequency, or scale-to-period conversions, as well as the cone % of influence. In this example, the helper function |helperPlotCoherence| % packages some useful commands for plotting the outputs of |wcoherence|. [wcoh,~,F,coi] = wcoherence(NIRSData(:,1),NIRSData(:,2),10,'numscales',16); helperPlotCoherence(wcoh,tm,F,coi,'Seconds','Hz'); %% % In the plot, you see a region of strong coherence throughout the data % collection period around 1 Hz. This results from the cardiac rhythms of % the two subjects. Additionally, you see regions of strong coherence % around 0.13 Hz. This represents coherent oscillations in the subjects' % brains induced by the task. If it is more natural to view the wavelet % coherence in terms of periods rather than frequencies, you can use the % 'dt' option and input the sampling interval. With the 'dt' option, % |wcoherence| provides scale-to-period conversions. [wcoh,~,P,coi] = wcoherence(NIRSData(:,1),NIRSData(:,2),seconds(0.1),... 'numscales',16); helperPlotCoherence(wcoh,tm,seconds(P),seconds(coi),... 'Time (secs)','Periods (Seconds)'); %% % Again, note the coherent oscillations corresponding to the subjects' % cardiac activity occurring throughout the recordings with a period of % approximately one second. The task-related activity is also apparent with % a period of approximately 8 seconds. Consult Cui, Bryant, & Reiss (2012) % for a more detailed wavelet analysis of this data. %% Conclusions % In this example you learned how to use wavelet coherence to look for % time-localized coherent oscillatory behavior in two time series. For % nonstationary signals, it is often more informative if you have a measure % of coherence that provides simultaneous time and frequency (period) % information. The relative phase information obtained from the wavelet % cross-spectrum can be informative when one time series directly affects % oscillations in the other. %% 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), pp. 2430-2437, 2012. % % Grinsted, A., Moore, J.C., and Jevrejeva, S. "Application of the cross % wavelet transform and wavelet coherence to geophysical time series", % Nonlin. Processes Geophys., 11, pp. 561-566, 2004. % % Maraun, D., Kurths, J. and Holschneider, M. "Nonstationary Gaussian % processes in wavelet domain: Synthesis, estimation and significance % testing", Phys. Rev. E 75, pp. 016707(1)-016707(14), 2007. % % Torrence, C. and Webster, P. "Interdecadal changes in the ESNO-Monsoon % System," J.Clim., 12, pp. 2679-2690, 1999. %% Appendix % The following helper function is used in this example. % % * <matlab:edit('helperPlotCoherence.m') helperPlotCoherence>