CWT operation command does not work as defined by Help Topics and Command Help

I have been wondering how to get correct CWT result.
At long last, I got near answer.
Please help me to get perfect answer.

ThreadSafe Static Function TSCwtMorlet(wSrc, wOut, vIdx, vM, WPR1, vScl, vCmplx)
 Wave wSrc, 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 /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=4 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, WPR1, vScl, vCmplx, vPrec)
 Wave wSrc, 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))
 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
  M_CWTI = r2polar(M_CWTI)
  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


Function/S MTCwtMorlet(sWave, vM, vCmplx, WPR1, SMP2, startScale, scaleStepSize, numScales, sW2)
 String sWave, sW2
 Variable vM, vCmplx, WPR1, SMP2, startScale, scaleStepSize, numScales
 
 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)
 Make /N=(vSize, numScales) $sOut
 Wave wOut=$sOut
 
 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)
     ThreadStart vMT, i, TSCwtMorlet(wSrc, wOut, vIdx, vM, WPR1, vScl, vCmplx)
    else
     ThreadStart vMT, i, TSCwtDefMorlet(wSrc, wOut, vIdx, WPR1, vScl, vCmplx, 0)
    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)
     TSCwtMorlet(wSrc, wOut, vIdx, vM, WPR1, vScl, vCmplx)
    else
     TSCwtDefMorlet(wSrc, wOut, vIdx, WPR1, vScl, vCmplx, 0)
    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)
 String sWave
 Variable vCmplx, WPR1=2*pi, SMP2, startScale=2, scaleStepSize=.5, numScales=12
 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"
 
 vCmplx = (vCmplx==1)
 SMP2 = SMP2 == 1 ? 1:4
 Variable vM = 0
 String sOut=MTCwtMorlet(sWave, vM, vCmplx, WPR1, SMP2, startScale, scaleStepSize, numScales, "")
 String sW2=StringFromList(1, sOut)
 do
  vM += 1
  MTCwtMorlet(sWave, vM, vCmplx, WPR1, 2, startScale, scaleStepSize, numScales, sW2)
 while (vM < 2)
end
It would be more helpful if you told us specifically what does not work for you.

A.G.
WaveMetrics, Inc.
Igor wrote:
It would be more helpful if you told us specifically what does not work for you.

A.G.
WaveMetrics, Inc.



Thank you for your attention.

CWT command gives me different results between fft method(/M=0) and direct integration method(/M=1).
Both result seems to be different with definition.
Direct integration method has little difference but it takes very long time.
FFT method answers relatively very fast but the result has large differece almost proportional to source wave size.
Because I am not familiar with English, I describe macro code.
CwtTest() makes three 2-dimentional waves.
The named wave as "*Def*" will be a result defined by Help Topics and Command Help.
The rest two are related with CWT command.
kanekaka wrote:

Because I am not familiar with English, I describe macro code.
CwtTest() makes three 2-dimentional waves.
The named wave as "*Def*" will be a result defined by Help Topics and Command Help.
The rest two are related with CWT command.


Thanks for the explanation. Now I can narrow down the specifics of your question. If I understand you correctly, it should be possible to illustrate the difference by executing 4 simple instructions:

1. Create a particular source data:
Make/n=1000 ddd=some function
2. Execute CWT with one set of flags
CWT [flags1] ddd
3. Store the output in another name:
Duplicate/O M_CWT ddd_flags1
4. Execute the CWT with a different set of flags
CWT[flags2] ddd

If you can send me such a set of instructions, I'd be able to investigate your report further. Also note that you can communicate this directly to support@wavemetrics.com.

A.G.
WaveMetrics, Inc.