How to interpret each row of the output matrix of CWT?
8 views (last 30 days)
Show older comments
Hi, I just started practicing CWT function included in the Wavelet Toolbox, and I encountered a very weird situation.
Probably I totally misunderstood the CWT function because of lack of mathematical background, so.... I'd like to ask you some help.
I got an output matrix W from a time series X (complex number) by using CWT function as below:
N = 1000; % number of sampling
dt = 1; % time interval: day
t = [0:N-1]*dt; % time vector: day
T = 60; % period: day
% sample time series with a single freq
X = 3*exp(1i*2*pi*1/T*t); % positive rotation on r=3 circle
% CWT
[W, period, coi] = cwt(X, days(1));
P = days(period);
% I found the P(48) is the closest period to T, and 48th row of abs(W) is
% also nearly equal to 3, the amplitude I assumed above.
j = 48;
Xj = W(j, :, 1); % pick the row (W(:,:,2) were almost zero in this example)
% insert the row into an empty matrix
W2 = zeros(size(W));
W2(j, :, 1) = Wj;
% inverse transform using W2
X2 = icwt(W2);
And the real parts of the results are like this:
I think the X2 from the inverse transform of a single row has a reduced amplitude because I ignored other "spreading" signals in W (like 47th and 49th rows). Is this correct?
And Xj includes a small boundary effect at both ends, but shows almost equal amplitudes to X that I intended.
So I thought that each j-th row of matrix W represented the oscillatory component of X at the j-th frequency, and that X3 = sum(W(:,:,1), 1) should be almost equal to X.
But of course, it was not; X4 = icwt(W) was the closest solution to X.
Here I found a very weird pattern:
The ratio Xj ./ X2 was always about 4.7273... (blue line in the third panel), and X3./X also converges to the same value (red).
I also tested other cases changing the length of data, or using more complicated time series, but I get the same ratio from Xj/X2 and X3/X.
What is this number come from? And why is this number constant?!
Is this like a scaling factor of 2/length(data) in FFT?
I really appreciate if anyone can explain the meaning of each row of W and scaling factor of inverse transform....
Thank you in advance.
0 Comments
Accepted Answer
Angelo Yeo
on 9 Oct 2024
Here is the code to reproduce your issue. (Next time, please put all the script to help others reproduce the issue)
N = 1000; % number of sampling
dt = 1; % time interval: day
t = [0:N-1]*dt; % time vector: day
T = 60; % period: day
% sample time series with a single freq
X = 3*exp(1i*2*pi*1/T*t); % positive rotation on r=3 circle
% CWT
[W, period, coi] = cwt(X, days(1));
P = days(period);
% I found the P(48) is the closest period to T, and 48th row of abs(W) is
% also nearly equal to 3, the amplitude I assumed above.
j = 48;
Xj = W(j, :, 1); % pick the row (W(:,:,2) were almost zero in this example)
% insert the row into an empty matrix
W2 = zeros(size(W));
W2(j, :, 1) = Xj;
% inverse transform using W2
X2 = icwt(W2);
X3 = sum(W(:,:,1), 1);
X4 = icwt(W);
% X2 = icwt(W);
%%
figure
subplot(3,1,1);
plot(real(X))
hold on
plot(real(Xj))
plot(real(X2))
legend(["X", "Xj", "X2"])
subplot(3,1,2);
plot(real(X));
hold on;
plot(real(X3));
plot(real(X4));
legend(["X", "X3", "X4"])
subplot(3,1,3);
plot(real(Xj./X2));
hold on;
plot(real(X3./X));
plot(real(X4./X));
legend(["Xj/X2", "X3/X", "X4/X"])
For the first question, yes, you ignored neighboring components of CWT results and reconstructed the signal. So, not all the wanted components of signal are reconstructed.
For your second question, when calculating the inverse Wavelet Transform in MATLAB, Morlet's single integral formula is used, and the final coefficient calculated as shown in the figure below, the reciprocal of 2*log(a0)*(1/cpsi) is 4.7273. (The reciprocal is calculated because you asked about the result of calculating it with the inverse result in the denominator and the original signal in the numerator.)
However, this coefficient is a value that can change depending on conditions. In other words, this scaling factor changes depending on how the wavelet is calculated. It is not a fixed value like the inverse Fourier transform.
You can probably find Morlet's single integral formula used in MATLAB's icwt function in many places, but you can also find an explanation of the coefficients of reconstruction in Mallat's "A wavelet tour of signal processing," which is the most famous wavelet textbook, as shown below.
0 Comments
More Answers (0)
See Also
Categories
Find more on Continuous Wavelet Transforms in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!