CWT operation command does not work as defined by Help Topics and Command Help
kanekaka
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
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
A.G.
WaveMetrics, Inc.
October 11, 2013 at 01:53 pm - Permalink
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.
October 11, 2013 at 03:50 pm - Permalink
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.
October 11, 2013 at 05:15 pm - Permalink