Frequency Scale in Wavelet Transform

Dear All, 

I am trying to apply continuous wavelet transform to sound data but the periods on y-axis are way off from those resulting from FFT. 

I tried to use CWT demo but the problem still existed.  

In the attachment, SigL wave is the data to be CWTed.  Signal is the original sound file. 

Any suggestions? 

Best,

IGOR file with sound data

Have you used the /FSCL flag?

I don't think that the mapping between the CWT and frequency analysis via the FFT is very rigorous so some caution is advised.

 

A.G.

The reason that I feel uncomfortable with this mapping is that IMO it is completely inappropriate.  The conventional notion of frequency applies to periodic basis functions such as sines and cosines.  The CWT wavelets do not represent pure frequencies so while I added support for the scaling factor to match its use in the literature, I don't think it is appropriate to display the vertical axis in units of e.g., Hz.  In reality, the vertical axis represents wavelet scale.  If exact frequency representation is important to you I recommend using the STFT or the WignerTransform.

AG

Dear Sorasak Danworaphong

 

Because frequency is reciprocal to scale parameter and y-axis of M_CWT is related to the scale, /SMP2=4 and display it by log-axis is suitable.

See the attachment file which was calculated from your file using following macro codes.

In order to decrease the file size some waves and objects are deleted.

 

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

 

Forum

Support

Gallery

Igor Pro 9

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More