bug in cwt function

I am using cwt function with 'morl' wavelet. Here is the syntax I am using:

[cwtwavex,freqsx]=cwt(x, scales, 'morl' , 'absglb', dt)

But it seems that there is a bug in the function, because when I choose an impulse function as the input (x), and plot the output in different scales, which should be essentially the wavelets shape, but the output on the last scale is mainly a morlet function modulated with another function!... Here is the code I wrote: clear all;

xlen=1000; Fs=2e4; % sampling frequency

dt=1/Fs; a0 = 2^(1/32); scales = 2*a0.^(0:6*32);

x=zeros(1,xlen+1);

x(int16(xlen/2))=1;

wavex=x; [cwtwavex,freqsx]=cwt(x, scales, 'morl' , 'absglb', dt);

abs_cwt=abs(cwtwavex);

%% figure(1); plot(x);

figure (2); plot(cwtwavex(1,:), 'r') ylabel('cwtx1');

figure (3); plot(cwtwavex(193,:), 'b') ylabel('cwtx193');

figure(4); mesh(abs_cwt);

But as you see, the shape of the wavelet output on the last scale (193) and also the mesh figure are very strange ( many fluctuations) which would stem from the wrong shape of the morlet wavelet at S(193) (last scale)!!

This problem does not exist when 'sym2' is used! Please help.

Also, it seems that there are lots of limitation when using morlet wavelet with cwt function: it doesn't take 'morl' easily with regular syntax. Please advise.

Thank you very much.

Best

Answer 1

[Wayne King]

Hi Ziba, First of all, I would encourage you to use the new CWT interface as opposed to the legacy interface you are using here. With the new CWT interface, you can use the 'amor' for the analytic Morlet.

As pertains to your current question, I'm not sure what you are expecting. Let's use a delta function so when we obtain the CWT, we essentially get back the wavelet at each scale (subject to a reversal).

Fs=2e4; % sampling frequency
dt=1/Fs; 
a0 = 2^(1/32); 
scales = 2*a0.^(0:6*32);
x = zeros(1e3,1);
x(500) = 1;
[cwtwavex,freqsx]=cwt(x, scales,'morl',dt);

Let's know check the frequency at the 193-th element of the scale vector

freqsx(193)

So that is 126 Hz. Now get the wavelet corresponding to that scale.

psi = cwtwavex(193,:);

Let's get the power spectrum, we expect the peak frequency to be around 127.

[Pxx,F] = periodgram(psi,[],2048,2e4);
[~,idx] = max(Pxx);
F(idx)

If we take a much smaller scale, we expect that the wavelets have a higher center frequency.

[Ziba Ebra]

Dear Wayne

Thank you very much for your answer. But I would like to plot wavelet function at different scales to make sure that they are correct. When you plot(psi) (in your example), it has not a correct shape for a scaled morlet!!, for it seems to be modulated with another function, and I believe that something should be wrong. Please help. Thank you very much.

Best Ziba

Answer 2

[Ziba Ebra]

Dear Wayne,

I examined the cwt function and my data again. "cwt" Matlab function has definitely bug and when applied on the data, many problems (which mostly is demonstrated as fluctuation ) occurs which causes deadly impact on the output!.... as said I would use morlet wavelet, with 6*32 scales (no of octave=6). You can test it yourself. It took me a good amount of time (2-3 months) before I finally realized that too much problem with my output is not from instrumentation or the signal itself, but the CWT MATLAB FUNCTION!!....I had to finally write a function by myself to do continuous wavelet analysis!!....

Hope that mathwork team correct this bug!!...

Best Ziba