Question for CWT operation

Share icon

Hallo.

I would like to analyze some time series with the continuous wavelet transform. As reference I use the journal article from Torrence and Compo (A Practical Guide to Wavelet Analysis).
The problem is that the results from CWT operation in Igor dose not match the results from the article and the provided software. Using the formulas I was unable to reproduce the results form Igor too.
Hence it would be helpful to know what the CWT operation exactly did and which formulas are used.

Thank you.
You are not really telling us what is "the provided software" nor what results you are trying to reproduce. As a matter of practicality I suggest that you pick up a simple analytic function, generate its time series and compute the transform. The result is going to be independent of ANY publication or external software.

You might also want to consult the relevant demo which you can find under File Menu->Example Experiments->Analysis->CWT Demo and the Time-Frequency tools (http://www.igorexchange.com/project/TFPlot), although the latter requires a simple change in one procedure file.

A.G.
WaveMetrics, Inc.

October 5, 2015 at 01:05 pm - Permalink

The software is available from the authors URL (http://paos.colorado.edu/research/wavelets/) in my case the Matlab version.

I have created the following time series:
Wave1
-0.282106
0.932004
1.02987
-0.0784919
-1.37405
-1.66569
-0.685259
0.665821
1.14539
0.312517

All four tables are calculated with the same parameters. Here I present only the first column (constant scaling and varying translation parameter) of the transformation.
The first table shows the transformation with the CWT operation only. The setting was: CWT/ENDM=0/OUT=1/Q/R2={2,0.01,300}/SMP2=4/WBI1=MorletC/WPR1=6 Wave1
The second table shows the results from Time-Frequency Toolkit.
The third shows the results from Torrence and Compo (Matlab) software.
In the fourth table I calculated the transformation manually with the formulas from the (Torrence and Compo) article.
Table.jpg (156.12 KB)

October 6, 2015 at 02:46 am - Permalink

peeda85 wrote:
The software is available from the authors URL (http://paos.colorado.edu/research/wavelets/) in my case the Matlab version.

I have created the following time series:
Wave1
-0.282106
0.932004
1.02987
-0.0784919
-1.37405
-1.66569
-0.685259
0.665821
1.14539
0.312517

All four tables are calculated with the same parameters. Here I present only the first column (constant scaling and varying translation parameter) of the transformation.
The first table shows the transformation with the CWT operation only. The setting was: CWT/ENDM=0/OUT=1/Q/R2={2,0.01,300}/SMP2=4/WBI1=MorletC/WPR1=6 Wave1
The second table shows the results from Time-Frequency Toolkit.
The third shows the results from Torrence and Compo (Matlab) software.
In the fourth table I calculated the transformation manually with the formulas from the (Torrence and Compo) article.


Please send a copy of this experiment to support@wavemetrics.com. It is difficult to investigate what is going on here without having access to the exact data and commands used.

A.G.
WaveMetrics, Inc.

October 6, 2015 at 12:54 pm - Permalink

Igor][quote=kanekaka]I also think that CWT operation is not correctly implemented.<br> <br> I've posted following article.<br> <br> <a href="http://www.igorexchange.com/node/5404[/quote">http://www.igorexchange.com/node/5404[/quote</a> wrote:

Our support database does not indicate that you followed up on that thread. Feel free to do so now.

A.G.
WaveMetrics, Inc.


Dear, A.G.

I want to say twice that problems is your matter.
At 1st, I have wanted to tell that CWT operation command includes defection. And I indicated ideas of correction method that will help to investigate what is wrong.
Next, your e-mail receiving server is so confused that he rejects your customers e-mail as judging UCE (Spam). You should tell your unskillful postmaster that he tarnishes companies brand. It is miserable as computer engineer company.


Dear, peeda85

Your data of Wave1 is too short to get meaningful CWT result.
You should better to acquire more precise or/and longer wave.

October 9, 2015 at 12:54 am - Permalink

kanekaka wrote:
Dear, peeda85

Your data of Wave1 is too short to get meaningful CWT result.
You should better to acquire more precise or/and longer wave.



Hallo kanekaka,

thank you for your experiment. Good to see what other people did.
I know that the time series is too short to make sense. Because I computed it also manually it was enough for a comparison.

October 8, 2015 at 04:27 am - Permalink

peeda85 wrote:
I know that the time series is too short to make sense. Because I computed it also manually it was enough for a comparison.


"Not make sense" means that any algorithm can not calculate 'true' answer but far approximation solution. So it is natural that different algorithm gave different result.

I modified my macro code in order to answer complex format.
MTCwtMorlet("Wave1", 2, 1, 6, 4, 2, .01, 300, 1, "") will give same result with your manual transformation. 

ThreadSafe Static Function TSCwtMorlet(wSrc, wOut, vIdx, vM, WPR1, vScl, vCmplx, OUT)
 Wave wSrc, wOut
 Variable vIdx, vM, WPR1, vScl, vCmplx, OUT
 
 Variable vSize = numpnts(wSrc), StX = pnt2x(wSrc, 0), DtX = deltax(wSrc)
 String sUnit = WaveUnits(wSrc, 0), WBI1 = StringFromList(vCmplx, "Morlet;MorletC")
 
 Make /O/N=(vSize,1) M_CWT
 Variable V_flag
 
 if (vM == 0)           // In FFT method, deltax(wSrc ) should be 1
  SetScale/P x 0, 1, sUnit, wSrc
 endif
 
 CWT /M=(vM)/WBI1={$WBI1}/WPR1=(WPR1)/R2={vScl, 1, 1 }/OUT=(OUT) wSrc
 
 if (vM == 0)           // FFT
  SetScale/P x StX, DtX, sUnit, wSrc
  FastOP M_CWT = (1/(2^ceil(ln(vSize)/ln(2)))) * M_CWT
 else               // Direct
  FastOP M_CWT = (sqrt(2*pi)/pi^.25) * M_CWT
 endif
 
 wOut[][vIdx] = M_CWT[p]
 KillWaves/Z M_CWT
 Wave/Z M_CWT = $""
 
 Return V_flag
end
 
 
ThreadSafe Static Function TSCwtDefMorlet(wSrc, wOut, vIdx, OUT, WPR1, vScl, vCmplx, vPrec)
 Wave wSrc, wOut
 Variable vIdx, OUT, WPR1, vScl, vCmplx, vPrec
 
 if (vPrec == 0)
  vPrec = 1e-4
 endif
 Variable vWid = ceil(sqrt(-2*ln(vPrec)))
 Variable vNum = ceil(vScl*vWid)
 
 String sTmp = NameOfWave(wSrc) + "MCWT" + num2str(round(vScl*100))
 Duplicate/O wSrc, $sTmp
 
 Variable vA = 1/(pi^.25)/sqrt(vScl), vB = -.5/vScl^2, vC = WPR1/vScl
 String sCoef = NameOfWave(wSrc) + "Coef" + num2str(round(vScl*100))
 if (vCmplx == 1)
  Make/O/C/N=(vNum * 2 + 1) $sCoef
  Wave/C wCoefC = $sCoef
  SetScale /P x (-vNum), 1, "", wCoefC
  wCoefC = vA * exp(vB*x^2)*cmplx(cos(vC*x), sin(vC*x))
  Redimension /C $sTmp
  Wave/C M_CWTI = $sTmp
  Convolve wCoefC, M_CWTI
  if (OUT == 4)
   M_CWTI = r2polar(M_CWTI)
  endif
  Redimension /R M_CWTI
  Wave/C/Z M_CWTI = $""
  Wave M_CWT = $sTmp
 else
  Make/O/N=(vNum * 2 + 1) $sCoef
  Wave wCoef = $sCoef
  SetScale /P x (-vNum), 1, "", wCoef
  wCoef = vA * exp(vB*x^2)*cos(vC*x)
  Wave M_CWT = $sTmp
  Convolve wCoef, M_CWT
 endif
 
 wOut[][vIdx] = M_CWT[p+vNum]
 
 KillWaves/Z wCoef, wCoefC, M_CWT
 Wave/Z M_CWT = $""
 Wave/Z wCoef = $""
 
 Return vNum
end
 
 
ThreadSafe Static Function TSCwtCMorlet(wSrc, wOut, vIdx, vM, WPR1, vScl, vCmplx)
 Wave wSrc
 Wave/C wOut
 Variable vIdx, vM, WPR1, vScl, vCmplx
 
 Variable vSize = numpnts(wSrc), StX = pnt2x(wSrc, 0), DtX = deltax(wSrc)
 String sUnit = WaveUnits(wSrc, 0), WBI1 = StringFromList(vCmplx, "Morlet;MorletC")
 
 Make /C/O/N=(vSize,1) M_CWT
 Variable V_flag
 
 if (vM == 0)           // In FFT method, deltax(wSrc ) should be 1
  SetScale/P x 0, 1, sUnit, wSrc
 endif
 
 CWT /M=(vM)/WBI1={$WBI1}/WPR1=(WPR1)/R2={vScl, 1, 1 }/OUT=1 wSrc
 
 if (vM == 0)           // FFT
  SetScale/P x StX, DtX, sUnit, wSrc
  FastOP/C M_CWT = (cmplx(1/(2^ceil(ln(vSize)/ln(2))),0)) * M_CWT
 else               // Direct
  FastOP/C M_CWT = (cmplx(sqrt(2*pi)/pi^.25,0)) * M_CWT
 endif
 
 wOut[][vIdx] = M_CWT[p]
 KillWaves/Z M_CWT
 Wave/C/Z M_CWT = $""
 
 Return V_flag
end
 
 
ThreadSafe Static Function TSCwtDefCMorlet(wSrc, wOut, vIdx, WPR1, vScl, vCmplx, vPrec)
 Wave wSrc
 Wave/C wOut
 Variable vIdx, WPR1, vScl, vCmplx, vPrec
 
 if (vPrec == 0)
  vPrec = 1e-4
 endif
 Variable vWid = ceil(sqrt(-2*ln(vPrec)))
 Variable vNum = ceil(vScl*vWid)
 
 String sTmp = NameOfWave(wSrc) + "MCWT" + num2str(round(vScl*100))
 Duplicate/O wSrc, $sTmp
 
 Variable vA = 1/(pi^.25)/sqrt(vScl), vB = -.5/vScl^2, vC = WPR1/vScl
 String sCoef = NameOfWave(wSrc) + "Coef" + num2str(round(vScl*100))
  Make/O/C/N=(vNum * 2 + 1) $sCoef
  Wave/C wCoefC = $sCoef
  SetScale /P x (-vNum), 1, "", wCoefC
 if (vCmplx == 1)
  wCoefC = vA * exp(vB*x^2)*cmplx(cos(vC*x), sin(vC*x))
 else
  wCoefC = vA * exp(vB*x^2)*cmplx(cos(vC*x), 0)
 endif
  Redimension /C $sTmp
  Wave/C M_CWTI = $sTmp
  Convolve wCoefC, M_CWTI
 
 wOut[][vIdx] = M_CWTI[p+vNum]
 
 KillWaves/Z M_CWTI, wCoefC
 Wave/C/Z M_CWT = $""
 Wave/C/Z wCoefC = $""
 
 Return vNum
end
 
 
Function/S MTCwtMorlet(sWave, vM, vCmplx, WPR1, SMP2, startScale, scaleStepSize, numScales, OUT, sW2)
 String sWave, sW2
 Variable vM, vCmplx, WPR1, SMP2, startScale, scaleStepSize, numScales, OUT
 
 Wave wSrc=$sWave
 Variable vSize = numpnts(wSrc)
 
 if (SMP2 == 2)
  Wave wW2=$sW2
  numScales = numpnts(wW2)
 else
  sW2 = UniqueName(sWave + "Scl", 1, 0)
  Make /N=(numScales) $sW2
  Wave wW2=$sW2
  if (SMP2 == 1)
   wW2 = startScale + scaleStepSize * p
  else
   scaleStepSize = 2 ^ scaleStepSize
   wW2 = startScale * scaleStepSize ^ p
  endif
 endif
 
 String sM = StringFromList(vM, "Fft;Direct;Def")
 String sOut = UniqueName(sWave + sM + "CWT", 1, 0)
 if (OUT == 1)
  Make /C/N=(vSize, numScales) $sOut
  Wave/C wcOut=$sOut
 else
  Make /N=(vSize, numScales) $sOut
  Wave wOut=$sOut
 endif
 
 printf "%s;CWT,", time()
 Variable vIdx = 0, vScl
// Only work version 6 or later.
Variable vPC = ThreadProcessorCount
Variable i, vTGS = 0, vTGW = 100, vMT = ThreadGroupCreate(vPC)
if (vPC>1)
  do
   for(i=0; i < vPC; i+= 1)
    vScl = wW2[vIdx]
    printf "%d:%.2f,", vIdx, vScl
    if (vM < 2)
     if (OUT == 1)
      ThreadStart vMT, i, TSCwtCMorlet(wSrc, wcOut, vIdx, vM, WPR1, vScl, vCmplx)
     else
      ThreadStart vMT, i, TSCwtMorlet(wSrc, wOut, vIdx, vM, WPR1, vScl, vCmplx, OUT)
     endif
    else
     if (OUT == 1)
      ThreadStart vMT, i, TSCwtDefCMorlet(wSrc, wcOut, vIdx, WPR1, vScl, vCmplx, 0)
     else
      ThreadStart vMT, i, TSCwtDefMorlet(wSrc, wOut, vIdx, OUT, WPR1, vScl, vCmplx, 0)
     endif
    endif
    vIdx += 1
    if (vIdx >= numScales)
     break
    endif
   endfor
   do
    vTGS = ThreadGroupWait(vMT, vTGW)
   while (vTGS != 0)
  while (vIdx < numScales)
  vTGS = ThreadGroupRelease(vMT)
else
// If version 5, above code should be removed.
  do
    vScl = wW2[vIdx]
    printf "%d:%.2f,", vIdx, vScl
    if (vM < 2)
     if (OUT == 1)
      TSCwtCMorlet(wSrc, wcOut, vIdx, vM, WPR1, vScl, vCmplx)
     else
      TSCwtMorlet(wSrc, wOut, vIdx, vM, WPR1, vScl, vCmplx, OUT)
     endif
    else
     if (OUT == 1)
      TSCwtDefCMorlet(wSrc, wcOut, vIdx, WPR1, vScl, vCmplx, 0)
     else
      TSCwtDefMorlet(wSrc, wOut, vIdx, OUT, WPR1, vScl, vCmplx, 0)
     endif
    endif
    vIdx += 1
  while (vIdx < numScales)
endif   // If version 5, this line should be removed.
 print time()
 Return sOut + ";" + sW2
end
 
 
Macro CwtTest(sWave, vCmplx, WPR1, SMP2, startScale, scaleStepSize, numScales, vOUT)
 String sWave
 Variable vCmplx, WPR1=2*pi, SMP2=2, startScale=2, scaleStepSize=.5, numScales=12, vOUT=3
 Prompt sWave, "source wave",  popup, WaveList("*", ";", "DIMS:1")
 Prompt vCmplx, "Select mother wavelet", popup, "MorletC;Morlet"
 Prompt WPR1, "Angular frequency OMEGA"
 Prompt SMP2, "Scaling method", popup, "Linear;Power of 2"
 Prompt startScale, "Start scale"
 Prompt scaleStepSize, "Scale step size"
 Prompt numScales, "Number of scales"
 Prompt vOUT, "Output format", popup, "complex;real part;magnitude"
 
 vCmplx = (vCmplx==1)
 SMP2 = SMP2 == 1 ? 1:4
 vOUT = 2^(vOUT-1)
 Variable vM = 0
 String sOut=MTCwtMorlet(sWave, vM, vCmplx, WPR1, SMP2, startScale, scaleStepSize, numScales, vOUT, "")
 String sW2=StringFromList(1, sOut)
 do
  vM += 1
  MTCwtMorlet(sWave, vM, vCmplx, WPR1, 2, startScale, scaleStepSize, numScales, vOUT, sW2)
 while (vM < 2)
end

October 9, 2015 at 12:57 am - Permalink