Question for CWT operation
peeda85
Mon, 10/05/2015 - 07:53 am
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 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
I've posted following article.
http://www.igorexchange.com/node/5404
October 5, 2015 at 07:26 pm - Permalink
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.
October 6, 2015 at 02:46 am - Permalink
I've posted following article.
http://www.igorexchange.com/node/5404[/quote]
Our support database does not indicate that you followed up on that thread. Feel free to do so now.
A.G.
WaveMetrics, Inc.
October 6, 2015 at 12:51 pm - Permalink
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
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
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
"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.
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