www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/CWTTimeFrequencyExample.m

    %% Time-Frequency Analysis with the Continuous Wavelet Transform 
%
% This example shows how to use the continuous wavelet transform (CWT) to
% analyze signals jointly in time and frequency. The example discusses the
% localization of transients where the CWT outperforms the short-time
% Fourier transform (STFT). The example also shows how to synthesize
% time-frequency localized signal approximations using the inverse CWT.
%
% The CWT is compared with the STFT in a number of the examples. You must
% have the Signal Processing Toolbox(TM) to run the |spectrogram|
% examples.

%% Time-Frequency Analysis of Modulated Signals
% Load a quadratic chirp signal and show a plot of its spectrogram. The
% signal's frequency begins at approximately 500 Hz at t = 0, decreases to
% 100 Hz at t=2, and increases back to 500 Hz at t=4. The sampling
% frequency is 1 kHz.

load quadchirp;
fs = 1000;
[S,F,T] = spectrogram(quadchirp,100,98,128,fs);
helperCWTTimeFreqPlot(S,T,F,'surf','STFT of Quadratic Chirp','Seconds','Hz')
%%
% Obtain a time-frequency plot of this signal using the CWT. 
[cfs,f] = cwt(quadchirp,fs,'WaveletParameters',[14,200]);
helperCWTTimeFreqPlot(cfs,tquad,f,'surf','CWT of Quadratic Chirp','Seconds','Hz')

%% 
% The CWT with the bump wavelet produces a time-frequency analysis very
% similar to the STFT.
%
% Frequency and amplitude modulation occur frequently in natural signals.
% Use the CWT to obtain a time-frequency analysis of an echolocation pulse
% emitted by a big brown bat (Eptesicus Fuscus). The sampling interval is 7
% microseconds. Use the bump wavelet with 32 voices per octave. Thanks to 
% Curtis Condon, Ken White, and Al Feng of the Beckman Center at the 
% University of Illinois for the bat data and permission to use it in this 
% example.
load batsignal
t = 0:DT:(numel(batsignal)*DT)-DT;
[cfs,f] = cwt(batsignal,'bump',1/DT,'VoicesPerOctave',32);
helperCWTTimeFreqPlot(cfs,t*1e6,f./1000,'surf','CWT of Bat Echolocation',...
    'Microseconds','kHz')
%%
% Obtain and plot the STFT of the bat data.
[S,F,T] = spectrogram(batsignal,50,48,128,1/DT);
helperCWTTimeFreqPlot(S,T.*1e6,F./1e3,'surf','Bat Echolocation (STFT)',...
    'Microseconds','kHz')
%%
% For both the simulated and natural modulated signals, the CWT provides
% results similar to the STFT.

%% Detection of Transients in Oscillations Using the CWT
% There are certain situations in time-frequency analysis where the CWT can
% provide a more informative time-frequency transform than the short-time
% Fourier transform. One such situation occurs when the signal is corrupted
% by transients. The appearance and disappearance of these transients often
% has physical significance. Therefore, it is important to be able to
% localize these transients in addition to characterizing oscillatory
% components in the signal. To simulate this, create a signal consisting of
% two sine waves with frequencies of 150 and 200 Hz. The sampling frequency
% is 1 kHz. The sine waves have disjoint time supports. The 150-Hz sine
% wave occurs between 100 and 300 milliseconds. The 200-Hz sine wave occurs
% from 700 milliseconds to 1 second. Additionally, there are two transients
% at 222 and 800 milliseconds. The signal is corrupted by noise.

rng default;
dt = 0.001;
t = 0:dt:1-dt;
addNoise = 0.025*randn(size(t));
x = cos(2*pi*150*t).*(t>=0.1 & t<0.3)+sin(2*pi*200*t).*(t>0.7);
x = x+addNoise;
x([222 800]) = x([222 800 ])+[-2 2];
figure;
plot(t.*1000,x);
xlabel('Milliseconds'); ylabel('Amplitude');
%%
% Zoom in on the two transients to see that they represent disturbances in
% the oscillations at 150 and 200 Hz.
subplot(2,1,1)
plot(t(184:264).*1000,x(184:264));
ylabel('Amplitude')
title('Transients')
axis tight;
subplot(2,1,2)
plot(t(760:840).*1000,x(760:840));
ylabel('Amplitude')
axis tight;
xlabel('Milliseconds')

%%
% Obtain and plot the CWT using the analytic Morlet wavelet.
figure;
[cfs,f] = cwt(x,1/dt,'amor');
contour(t*1000,f,abs(cfs))
grid on
c = colorbar;
xlabel('Milliseconds')
ylabel('Frequency')
c.Label.String = 'Magnitude';
%%
% The analytic Morlet wavelet exhibits poorer frequency localization than
% the bump wavelet, but superior time localization. This makes the Morlet
% wavelet a better choice for transient localization. Plot the
% magnitude-squared fine scale coefficients to demonstrate the localization
% of the transients.

figure;
wt = cwt(x,1/dt,'amor');
plot(t.*1000,abs(wt(1,:)).^2);
xlabel('Milliseconds'); ylabel('Magnitude');
grid on; 
title('Analytic Morlet Wavelet -- Fine Scale Coefficients');
hold off
%%
% The wavelet shrinks to enable time localization of the transients with a
% high degree of accuracy while stretching to permit frequency localization
% of the oscillations at 150 and 200 Hz.

%%
% The STFT can only localize the transients to the width of the window. You
% have to plot the STFT in decibels (dB) to be able to visualize the
% transients.

[S,F,T] = spectrogram(x,50,48,128,1000);
surf(T.*1000,F,20*log10(abs(S)),'edgecolor','none'); view(0,90);
axis tight;
shading interp;
colorbar
xlabel('Time'), ylabel('Hz');
title('STFT')
%%
% The transients appear in the STFT only as a broadband increase
% in power. Compare short-time power estimates obtained from the STFT
% before (centered at 183 msec) and after (centered at 223 msec) the
% appearance of the first transient.

figure;
plot(F,20*log10(abs(S(:,80))));
hold on;
plot(F,20*log10(abs(S(:,100))),'r');
legend('T = 183 msec','T = 223 msec')
xlabel('Hz');
ylabel('dB');
hold off;
%%
% The STFT does not provide the ability to localize transients to the same
% degree as the CWT.
%% Removing A Time-Localized Frequency Component Using the Inverse CWT
% Create a signal consisting of exponentially weighted sine waves. There
% are two 25-Hz components -- one centered at 0.2 seconds and one centered
% at 0.5 seconds. There are two 70-Hz components -- one centered at 0.2 and
% one centered at 0.8 seconds. The first 25-Hz and 70-Hz components
% co-occur in time.

t = 0:1/2000:1-1/2000;
dt = 1/2000;
x1 = sin(50*pi*t).*exp(-50*pi*(t-0.2).^2);
x2 = sin(50*pi*t).*exp(-100*pi*(t-0.5).^2);
x3 = 2*cos(140*pi*t).*exp(-50*pi*(t-0.2).^2);
x4 = 2*sin(140*pi*t).*exp(-80*pi*(t-0.8).^2);
x = x1+x2+x3+x4;
figure;
plot(t,x)
title('Superimposed Signal')

%%
% Obtain and display the CWT.

cwt(x,2000);
title('Analytic CWT using Default Morse Wavelet');

%%
% Remove the 25 Hz component which occurs from approximately 0.07 to 0.3
% seconds by zeroing out the CWT coefficients. Use the inverse CWT
% (<matlab:doc('icwt') |icwt|>) to reconstruct an approximation to the
% signal.

[cfs,f] = cwt(x,2000);
T1 = .07;  T2 = .33;
F1 = 19;   F2 = 34;
cfs(f > F1 & f < F2, t> T1 & t < T2) = 0;
xrec = icwt(cfs);

%%
% Display the CWT of the reconstructed signal. The initial 25-Hz component
% is removed.

cwt(xrec,2000)

%%
% Plot the original signal and the reconstruction.

subplot(2,1,1);
plot(t,x);
title('Original Signal');
subplot(2,1,2);
plot(t,xrec)
title('Signal with first 25-Hz component removed');
%%
% Finally, compare the reconstructed signal with the original signal
% without the 25-Hz component centered at 0.2 seconds.
y = x2+x3+x4;
figure;
plot(t,xrec)
hold on
plot(t,y,'r--')
legend('Inverse CWT approximation','Original signal without 25-Hz');
hold off

%% Determining Exact Frequency Through the Analytic CWT
% When you obtain the wavelet transform of a sine wave using an analytic
% wavelet, the analytic CWT coefficients actually encode the frequency.
%
% To illustrate this, consider an otoacoustic emission obtained from a
% human ear. Otoacoustic emissions (OAEs) are emitted by the cochlea (inner
% ear) and their presence are indicative of normal hearing. Load and plot
% the OAE data. The data are sampled at 20 kHz.

load dpoae
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.

cwt(dpoaets,'bump',20000);
xlabel('Milliseconds');
%%
% 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.

[dpoaeCWT,f] = cwt(dpoaets,'bump',20000);
[~,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('Milliseconds'), ylabel('Magnitude');
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.

%Scidx = find(f >= 1150 & f <= 1350);
%emissionCWT = dpoaeCWT;
%tmpcfs = zeros(size(emissionCWT));
%tmpcfs(Scidx,:) = emissionCWT(Scidx,:);
%emissionCWT = tmpcfs;

xrec = icwt(dpoaeCWT,f,[1150 1350]);
figure
plot(t.*1000,dpoaets)
hold on;
plot(t.*1000,xrec,'r')
AX = gca;
ylim = AX.YLim;
plot([25 25],ylim,'k')
plot([175 175],ylim,'k')
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));
%% Conclusions and Further Reading
% In this example you learned how to use the CWT to obtain a time-frequency
% analysis of a 1-D signal using an analytic wavelet with |cwt|. You saw 
% examples of signals where the CWT provides similar results to the STFT 
% and an example where the CWT can provide more interpretable results than 
% the STFT. Finally, you learned how to reconstruct time-scale (frequency) 
% localized approximations to a signal using |icwt|.
%% 
% For more information on the CWT see the Wavelet Toolbox documentation.
%
% Reference: Mallat, S. "A Wavelet Tour of Signal Processing: The Sparse
% Way", Academic Press, 2009.
%% Appendix
% The following helper functions are used in this example.
%
% * <matlab:edit('helperCWTTimeFreqPlot.m') helperCWTTimeFreqPlot>
% * <matlab:edit('helperCWTTimeFreqVector.m') helperCWTTimeFreqVector>