Frequency Scale in Wavelet Transform

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) && (DtX != 1))           // In FFT method, deltax(wSrc) should be 1.
  SetScale/P x 0, 1, sUnit, wSrc        // This code does not executed here in fact.
 endif
 
 CWT /M=(vM)/WBI1={$WBI1}/WPR1=(WPR1)/R2={vScl, 1, 1 }/OUT=(OUT) wSrc
 
 if (vM == 0)               // FFT
  FastOP M_CWT = (1/(2^ceil(ln(vSize)/ln(2)))) * M_CWT  // M_CWT with FFT method will be large almost proportional to numpnts.
 else                       // Direct
  FastOP M_CWT = (sqrt(2*pi)/pi^.25) * M_CWT        // M_CWT with direct method should be devided with pi^(1/4) instead of sqrt(2 pi).
 endif
 
 wOut[][vIdx] = M_CWT[p]
 KillWaves/Z M_CWT
 Wave/Z M_CWT = $""
 
 Return V_flag
end
 
// CWT is defined as the convolution formula at each scale.
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
 
 // Wavelet with scale of vScl will be prepared as a wave of $sCoef.
 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
 
 CWT /M=(vM)/WBI1={$WBI1}/WPR1=(WPR1)/R2={vScl, 1, 1 }/OUT=1 wSrc
 
 if (vM == 0)               // FFT
  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), StX = pnt2x(wSrc, 0), DtX = deltax(wSrc)
 String sUnit = WaveUnits(wSrc, 0), sFW = ""
 
 Variable vF = WPR1/(2*pi), vL = 1 / DtX
 vL *= (vF+sqrt(vF^2+.5/pi^2))/2
 
 if (SMP2 == 2)
  Wave wW2=$sW2
  numScales = numpnts(wW2)
 else
  sW2 = UniqueName(sWave + "Scl", 1, 0)
  Make /N=(numScales + 1) $sW2
  Wave wW2=$sW2
  sFW = UniqueName(sWave + "YFreq", 1, 0)
  Make /N=(numScales + 1) $sFW      // Y wave for image plot should have DimSize(M_CWT, 1) + 1 points.
  Wave wFW=$sFW
  if (SMP2 == 1)
   wW2 = startScale + scaleStepSize * p
  else
   scaleStepSize = 2 ^ scaleStepSize
   wW2 = startScale * scaleStepSize ^ p
  endif
  if (SMP2 >= 3)
   wFW = wW2 / sqrt(scaleStepSize)
  else
   wFW = wW2 - .5*scaleStepSize
  endif
  wFW = vL / wFW                            // Fourier frequency is reciprocal number of Fourier wavelength lambda.
  if (!cmpstr(sUnit, "s"))
   SetScale d 0,0,"Hz", wFW
  endif
  Redimension /N=(numScales) $sW2
 endif
 
 if ((vM == 0) && (DtX != 1))           // In FFT method, deltax(wSrc) should be 1
  SetScale/P x 0, 1, sUnit, wSrc
 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
 SetScale/P x StX, DtX, sUnit, $sOut
 
 printf "%s;%sCWT,", time(), sM
 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.
 
 if ((vM == 0) && (DtX != 1))           // When FFT method and deltax(wSrc) was not 1, x scaling should be restored.
  SetScale/P x StX, DtX, sUnit, wSrc
 endif
 
 print time()
 Return sOut + ";" + sW2 + ";" + sFW
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