www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/synchrosqueezingExample.m
%% Time-Frequency Reassignment and Mode Extraction with Synchrosqueezing % This example shows how to use wavelet synchrosqueezing to obtain a higher % resolution time-frequency analysis. The example also shows how to % extract and reconstruct oscillatory modes in a signal. % % In many practical applications across a wide range of disciplines, % signals occur which consist of a number of oscillatory components, or % modes. These components often exhibit slow variations in amplitude and % smooth changes in frequency over time. Signals consisting of one or more % such components are called amplitude and frequency modulated (AM-FM). % Individual AM-FM components of signals are also referred to as intrinsic % modes, or intrinsic mode functions (IMF). % % Wavelet synchrosqueezing is a method for analyzing signals consisting of % a sum of well-separated AM-FM components, or IMFs. With synchrosqueezing, % you can sharpen the time-frequency analysis of a signal as well as % reconstruct individual oscillatory modes for isolated analysis. %% Sharpen Time-Frequency Analysis % Synchrosqueezing can compensate for the spreading in time and frequency % caused by linear transforms like the short-time Fourier and continuous % wavelet transforms (CWT). In the CWT, the wavelet acts like a measuring % device for the input signal. Accordingly, any time frequency analysis % depends not only on the intrinsic time-frequency properties of the % signal, but also on the properties of the wavelet. % % To demonstrate this, first obtain and plot the CWT of a quadratic chirp % signal. 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; [cfs,f] = cwt(quadchirp,'bump',fs); helperCWTTimeFreqPlot(cfs,tquad,f,'surf','CWT of Quadratic Chirp','Seconds','Hz') %% % Note that the energy of the quadratic chirp is smeared in the % time-frequency plane by the time-frequency concentration of the wavelet. % For example, if you focus on the time-frequency concentration of the CWT % magnitudes near 100 Hz, you see that it is narrower than that observed % near 500 Hz. This is not an intrinsic property of the chirp. It is % an artifact of the measuring device (the CWT). Compare the % time-frequency analysis of the same signal obtained with the % synchrosqueezed transform. wsst(quadchirp,1000,'bump') %% % The synchrosqueezed transform uses the phase information in the CWT to % sharpen the time-frequency analysis of the chirp. %% Reconstruct Signal from Synchrosqueezed Transform % You can reconstruct a signal from the synchrosqueezed transform. % This is an advantage synchrosqueezing has over other time-frequency % reassignment techniques. The transform does not provide a perfect % inversion, but the reconstruction is stable and the results are typically % quite good. To demonstrate this, load, plot, and play a recording of a % female speaker saying "I saw the sheep". load wavsheep; plot(tsh,sheep) title(' "I saw the sheep." '); xlabel('Time (secs)'); ylabel('Amplitude'); grid on; hplayer = audioplayer(sheep,fs); play(hplayer); %% % Plot the synchrosqueezed transform of the speech sample. [sst,f] = wsst(sheep,fs,'bump'); contour(tsh,f,abs(sst)); title('Wavelet Synchrosqueezed Transform'); xlabel('Time (secs)'); ylabel('Hz'); grid on; %% % Reconstruct the signal from the synchrosqueezed transform and compare the % reconstruction to the original waveform. xrec = iwsst(sst,'bump'); plot(tsh,sheep) hold on; plot(tsh,sheep,'r'); xlabel('Time (secs)'); ylabel('Amplitude'); grid on; title('Reconstruction From Synchrosqueezed Transform'); legend('Original Waveform','Synchrosqueezed Reconstruction'); sprintf('The maximum sample-by-sample difference is %1.3f', ... norm(abs(xrec-sheep),Inf)) hold off; %% % Play the reconstructed signal and compare with the original. hplayerRecon = audioplayer(xrec,fs); play(hplayerRecon) %% % The ability to reconstruct from the synchrosqueezed transform enables you % to extract signal components from localized regions of the time-frequency % plane. The next section demonstrates this idea. %% Identify Ridges and Reconstruct Modes % A time-frequency ridge is defined by the local maxima of a time-frequency % transform. Because the synchrosqueezed transform concentrates oscillatory % modes in a narrow region of the time-frequency plane and is invertible, % you can reconstruct individual modes by: % % # Identifying ridges in the magnitudes of the synchrosqueezed transform % # Reconstructing along the ridge % % This allows you to isolate and analyze modes which may be difficult or % impossible to extract with conventional bandpass filtering. % % To illustrate this, consider an echolocation pulse emitted by a big brown % bat (Eptesicus Fuscus). The sampling interval is 7 microseconds. 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 the data and plot the synchrosqueezed transform. load batsignal; time = 0:DT:(numel(batsignal)*DT)-DT; [sst,f] = wsst(batsignal,1/DT); contour(time.*1e6,f./1000,abs(sst)); grid on; xlabel('microseconds'); ylabel('kHz'); title('Bat Echolocation Pulse'); %% % Note that there are two modulated modes that trace out curved paths in % the time-frequency plane. Attempting to separate these components with % conventional bandpass filtering does not work because the filter would % need to operate in a time-varying manner. For example, a conventional % filter with a passband of 18 to 40 kHz would capture the energy in the % earliest-occurring pulse, but would also capture energy from the % later pulse. % % Synchrosqueezing can separate these components by filtering and % reconstructing the synchrosqueezed transform in a time-varying manner. % % First, extract the two highest-energy ridges from the synchrosqueezed % transform. [fridge,iridge] = wsstridge(sst,5,f,'NumRidges',2); hold on; plot(time.*1e6,fridge./1e3,'k--','linewidth',2); hold off; %% % The ridges follow the time-varying nature of the modulated pulses. % Reconstruct the signal modes by inverting the synchrosqueezed transform. xrec = iwsst(sst,iridge); subplot(2,1,1) plot(time.*1e6,batsignal); hold on; plot(time.*1e6,xrec(:,1),'r'); grid on; ylabel('Amplitude'); title('Bat Echolocation Pulse with Reconstructed Modes'); legend('Original Signal','Mode 1','Location','SouthEast'); subplot(2,1,2); plot(time.*1e6,batsignal); hold on; grid on; plot(time.*1e6,xrec(:,2),'r'); xlabel('microseconds'); ylabel('Amplitude'); legend('Original Signal','Mode 2','Location','SouthEast'); %% % You see that synchrosqueezing has extracted the two modes. If you sum the % two dominant modes at each point in time, you essentially recover the % entire echolocation pulse. In this application, synchrosqueezing allows % you to isolate signal components where a traditional bandpass filter % would fail. %% Penalty Term in Ridge Extraction % The previous mode extraction example used a penalty term in the ridge % extraction without explanation. % % When you extract multiple ridges, or you have a single modulated % component in additive noise, it is important to use a penalty term in the % ridge extraction. The penalty term serves to prevent jumps in frequency % as the region of highest energy in the time-frequency plane moves. % % To demonstrate this, consider a two-component signal consisting of an % amplitude and frequency-modulated signal plus a sine wave. The signal is % sampled at 2000 Hz. The sine wave frequency is 18 Hz. The AM-FM signal is % defined by: % % $(2+\cos(4\pi t))\sin(2\pi 231t+90\sin(3\pi t))$ % % Load the signal and obtain the synchrosqueezed transform. load multicompsig; sig = sig1+sig2; [sst,f] = wsst(sig,sampfreq); figure; contour(t,f,abs(sst)); grid on; title('Synchrosqueezed Transform of AM-FM Signal Plus Sine Wave'); xlabel('Time (secs)'); ylabel('Hz'); %% % First attempt to extract two ridges from the synchrosqueezed transform % without using a penalty. [fridge,iridge] = wsstridge(sst,f,'NumRidges',2); hold on; plot(t,fridge,'k--','linewidth',2); %% % You see that the ridge jumps between the AM-FM signal and the sine wave % as the region of highest energy in the time-frequency plane changes % between the two signals. Add a penalty term of 5 to the ridge extraction. % In this case, jumps in frequency are penalized by a factor of 5 times the % distance between the frequencies in terms of bins (not actual frequency % in hertz). [fridge,iridge] = wsstridge(sst,5,f,'NumRidges',2); figure; contour(t,f,abs(sst)); grid on; title('Synchrosqueezed Transform of AM-FM Signal Plus Sine Wave'); xlabel('Time (secs)'); ylabel('Hz'); hold on; plot(t,fridge,'k--','linewidth',2); hold off; %% % With the penalty term, the two modes of oscillation are isolated in two % separate ridges. Reconstruct the modes along the time-frequency % ridges from the synchrosqueezed transform. xrec = iwsst(sst,iridge); %% % Plot the reconstructed modes along with the original signals for % comparison. subplot(2,2,1) plot(t,xrec(:,1)); grid on; ylabel('Amplitude'); title('Reconstructed Mode'); ylim([-1.5 1.5]); subplot(2,2,2); plot(t,sig2); grid on; ylabel('Amplitude'); title('Original Component'); ylim([-1.5 1.5]); subplot(2,2,3); plot(t,xrec(:,2)); grid on; xlabel('Time (secs)'); ylabel('Amplitude'); title('Reconstructed Mode'); ylim([-1.5 1.5]); subplot(2,2,4); plot(t,sig1); grid on; xlabel('Time (secs)'); ylabel('Amplitude'); title('Original Component'); ylim([-1.5 1.5]); %% Conclusions % In this example, you learned how to use wavelet synchrosqueezing to % obtain a higher-resolution time-frequency analysis. % % You also learned how to identify maxima ridges in the synchrosqueezed % transform and reconstruct the time waveforms corresponding to those % modes. Mode extraction from the synchrosqueezing transform can help you % isolate signal components, which are difficult or impossible to isolate % with conventional bandpass filtering. %% References % % Daubechies, I., Lu, J., and Wu, H-T. "Synchrosqueezed wavelet transforms: % an empirical mode decomposition-like tool." Appl. Comput. Harmonic Anal., % 30(2):243-261, 2011. % % Thakur, G., Brevdo, E., Fuckar, N.S. and Wu H-T. "The synchrosqueezing % algorithm for time-varying spectral analysis: robustness properties and % new paleoclimate applications." Signal Processing, 93(5), 1079-1094, % 2013. % % Meignen, S., Oberlin, T., and McLaughlin, S. "A new algorithm for % multicomponent signals analysis based on synchrosqueezing: with an % application to signal sampling and denoising." IEEE Transactions on % Signal Processing, vol. 60, no. 12, pp. 5787-5798, 2012. %% Appendix % The following helper function is used in this example. % % * <matlab:edit('helperScaleVector.m') helperScaleVector>