Documentation |
On this page… |
---|
CWT of Sum of Disjoint Sinusoids Approximate Scale-Frequency Conversions |
To implement a DFT-based continuous wavelet analysis in the MATLAB^{®} command window, use cwtft and icwtft.
For the mathematical basis of the DFT-based continuous wavelet analysis and synthesis see:
The signal is a sum of two disjoint sinusoids. The sampling frequency is 1023 Hz. The total signal duration is 1 second. The frequencies of the two sine waves are 4 and 8 Hz. The 4-Hz sine wave has support on the initial 1/2 second of the 1-second interval. The 8-Hz sine wave has support on the final 1/2 second.
N = 1024; t = linspace(0,1,N); dt =1/(N-1); Y = sin(8*pi*t).*(t<=0.5) + sin(16*pi*t).*(t>0.5);
Obtain the continuous wavelet transform (CWT) using the default analytic Morlet wavelet, and plot the results.
sig = {Y,dt};
cwtS1 = cwtft(sig,'plot');
The figure shows a plot of the original signal. The CWT moduli, real and imaginary parts of the CWT coefficients, and the CWT coefficient arguments (phase angles) also appear as plots.
You can display the reconstructed signal by enabling the radio button at the bottom-left corner of the plot. Enabling the radio button superimposes the reconstructed signal on the original signal in the top-left corner of the figure. The relative maximum and quadratic (L2 norm) errors appear under the plot.
You can customize your continuous wavelet analysis by providing additional inputs to cwtft. In the following example, specify the analyzing wavelet as the Paul wavelet of order 8. Specify the initial scale, the spacing between scales, and the number of scales. By default, the scale vector is logarithmic to the base 2.
% smallest scale, spacing between scales, number of scales dt = 1/1023; s0 = 2*dt; ds = 0.5; NbSc = 20; % scale vector is % scales = s0*2.^((0:NbSc-1)*ds); wname = 'paul'; SIG = {Y,dt}; % Create SCA input as cell array SCA = {s0,ds,NbSc}; % Specify wavelet and parameters WAV = {wname,8}; % Compute and plot the CWT cwtS2 = cwtft(SIG,'scales',SCA,'wavelet',WAV,'plot');
The oscillations at 4 and 8 Hz are clearly visible as alternating positive and negative real and imaginary parts. The 4-Hz oscillation occurs at a longer scale than the 8-Hz oscillation. In the plot of the CWT moduli, you see the transition from the 4-Hz (longer scale) sinusoid to the 8-Hz sinusoid (shorter scale) around 0.5 seconds.
There is not a direct correspondence between Fourier wavelength and scale. However, you can find conversion factors for select wavelets that yield an approximate scale-frequency correspondence. You can find these factors for wavelets supported by cwtft listed on the reference page.
This example shows you how to change the scale axis to an approximate frequency axis for analysis. Use the sum of disjoint sinusoids as the input signal. Set the initial scale to 6*dt, the scale increment to 0.15, and the number of scales to 50. Use the Fourier factor for the Morlet wavelet to convert the scale vector to an approximate frequency vector in hertz. Plot the result.
figure; s0 = 6*dt; ds = 0.15; NbSc = 50; wname = 'morl'; SCA = {s0,ds,NbSc}; cwtsig = cwtft({Y,dt},'scales',SCA,'wavelet',wname); MorletFourierFactor = 4*pi/(6+sqrt(2+6^2)); Scales = cwtsig.scales.*MorletFourierFactor; Freq = 1./Scales; imagesc(t,[],abs(cwtsig.cfs)); indices = get(gca,'ytick'); set(gca,'yticklabel',Freq(indices)); xlabel('Time'); ylabel('Hz'); title('Time-Frequency Analysis with CWT');
You can see the signal contains significant energy at approximately 4 Hz over the first 1/2 second. In the final 1/2 second interval, the predominant signal energy transitions higher in frequency to approximately 8 Hz.
Repeat the above example using the Paul analyzing wavelet with order, m, equal to 8. Use a contour plot of the real part of the CWT to visualize the sine waves at 4 and 8-Hz. The real part exhibits oscillations in the sign of the wavelet coefficients at those frequencies.
s0 = 6*dt; ds = 0.15; NbSc = 50; m = 8; % scale vector is % scales = s0*2.^((0:NbSc-1)*ds); wname = 'paul'; SIG = {Y,dt}; % Create SCA input as cell array SCA = {s0,ds,NbSc}; % Specify wavelet and parameters WAV = {wname,m}; cwtPaul = cwtft(SIG,'scales',SCA,'wavelet',WAV); scales = cwtPaul.scales; PaulFourierFactor = 4*pi/(2*m+1); Freq = 1./(PaulFourierFactor.*scales); contour(t,Freq,real(cwtPaul.cfs)); xlabel('Time'); ylabel('Hz'); colorbar; title('Real Part of CWT using Paul Wavelet (m=8)'); axis([0 1 1 15]); grid on;
You can use the critically sampled (decimated) and oversampled (nondecimated) discrete wavelet transforms (DWT) to achieve perfect reconstruction of the input signal from the wavelet coefficients. To obtain a time and scale-dependent approximation to a signal, you can use a possibly-modified subset of the decimated or undecimated DWT coefficients.
The inversion of the CWT is not as straightforward. The simplest CWT inversion utilizes the single integral formula due to Morlet, which employs a Dirac delta function as the synthesizing wavelet. See Inverse Continuous Wavelet Transform for a brief mathematical motivation. icwtft and icwtlin both implement the single integral formula. Because of necessary approximations in the implementation of the single integral inverse CWT, you cannot expect to obtain perfect reconstruction. However, you can use the inverse CWT to obtain useful position and scale-dependent approximations to the input signal.
Implement the inverse CWT with logarithmically-spaced scales.
N = 1024; t = linspace(0,1,N); dt =1/(N-1); Y = sin(8*pi*t).*(t<=0.5) + sin(16*pi*t).*(t>0.5); dt = 1/1023; s0 = 2*dt; ds = 0.5; NbSc = 20; % scale vector is % scales = s0*2.^((0:NbSc-1)*ds); wname = 'paul'; SIG = {Y,dt}; % Create SCA input as cell array SCA = {s0,ds,NbSc}; % Specify wavelet and parameters WAV = {wname,8}; cwtS2 = cwtft(SIG,'scales',SCA,'wavelet',WAV); YR1 = icwtft(cwtS2,'plot','signal',SIG); norm(Y-YR1,2)
Enable the radio button in the left corner of the figure to plot the reconstructed signal.
Obtain the CWT of a noisy Doppler (frequency-modulated) signal using the analytic Morlet wavelet. Reconstruct an approximation by selecting a subset of the CWT coefficients. By eliminating the smallest scales, you obtain a lowpass approximation to the signal. The lowpass approximation produces a smooth approximation to the lower-frequency features of the noisy Doppler signal. The high-frequency (small scale) features at the beginning of the signal are lost.
% Define the signal load noisdopp; Y = noisdopp; N = length(Y); % Define parameters before analysis dt = 1; s0 = 2*dt; ds = 0.4875; NbSc = 20; wname = 'morl'; SIG = {Y,dt}; SCA = {s0,ds,NbSc}; WAV = {wname,[]}; % Compute CWT analysis cwtS4 = cwtft(SIG,'scales',SCA,'wavelet',WAV); % Thresholding step building the new structure cwtS5 = cwtS4; newCFS = zeros(size(cwtS4.cfs)); newCFS(11:end,:) = cwtS4.cfs(11:end,:); cwtS5.cfs = newCFS; % Reconstruction from the modified structure YRDen = icwtft(cwtS5,'signal',SIG); plot(Y,'k-.'); hold on; plot(YRDen,'r','linewidth',3); axis tight; legend('Original Signal', 'Selective inverse CWT'); title('Signal approximation based on a subset of CWT coefficients');