## Analysis of jumps in a two-level system

Average rating

Inspired by a recent discussion, here's a little snippet to analyze jumps in a two-level system, represented by (possibly very noisy) data that are distributed around a value A (for state A) and, from time to time, a value B. More documentation is in the snippet...

```// function AnalyzeTLS(inWave, decimation, minLifeTime, verbosity)
//	analyzes jumps in a two-level system
//
// 	Given an _inWave_ with the following characteristics:
//		- noisy data with two levels A (ground state) and B > A (excited state) (see #1 in code)
//		- both levels have approximately gaussian distribution (see assumption #2 in code)
//		- the ground state is occupied more often than the excited state (see #3 in code)
//	TwoStateAnalyze finds:
//		- where the jumps between A and B occur (output: EventStart)
//		- how long B is occupied in each case (output: LifeTime)
//		- how high the B level is above the mean A level in each case (output: EventLevel)
//
// 	Details and input parameters of AnalyzeTLS:
//		-  first decimates the data to reduce noise
//			- set _decimation_ factor to 1 for no decimation
//		- then analyzes the histogram to find mean levels for A and B
//		- uses these mean levels to find (usually too many) candidate events via FindLevels
//		- uses a lifetime criterion to reject spurious events (too short)
//			- set _minLifeTime_ to the desired value (e.g., 4)
//		- prints results to history (_verbosity_ > 0)
//		- displays histogram and inWave, overlaid with colored events (_verbosity_ > 1)
//
//	Original idea and first version by Xia Liu (Canada)
//	Rewritten by W. Harneit (Berlin, Germany) --- 2010/03/25
//
wave inWave			// wave to analyze
variable decimation	// factor for downscaling original data (> 1)
variable verbosity		// v = 0: be quiet, v = 1: print results but no graphs
// v > 1: also show histogram and colored data
DFREF saveDF = GetDataFolderDFR()
NewDataFolder/O/S root:Packages
NewDataFolder/O/S TLSAnalysis
DFREF pkgDF = GetDataFolderDFR()
// step 1: downsample data (eliminates noise)
string sampWaveName = NameOfWave(inWave)+"_samp"
Duplicate/O inWave \$sampWaveName
wave sampWave = \$sampWaveName
Resample/DOWN=(decimation) sampWave
// step 2: create histogram
string histWaveName = NameOfWave(inWave)+"_hist"
Make/N=100/O \$histWaveName
wave histWave = \$histWaveName
Histogram/B=1 sampWave, histWave
if( verbosity > 1 )
DoWindow/K histGraph			// kill histogram graph if it exists
Display/N=histGraph histWave
ModifyGraph log(left)=1
endif
// step 3: analyze histogram
WaveStats/M=1/Q histWave
variable levelA = V_maxLoc, valA = V_max
variable wid = levelA - leftx(histWave) 	// assume left (#1) peak is at edge (#2), and absolute maximum (#3)
WaveStats/M=1/Q/R=(levelA + wid, inf) histWave	// 2nd peak should be in this range
variable levelB = V_maxLoc, valB = V_max
if( verbosity > 0 )
printf "mean level(A) = %g (%g), mean level(B) = %g (%g)\r", levelA, valA, levelB, valB
endif
// step 4: curve fit histogram	// MAY BE UNNECESSARY
Duplicate/O histWave, \$(histWaveName+"_wt")
wave histWaveWeight = \$(histWaveName+"_wt")
histWaveWeight = sqrt(histWaveWeight)
Make/D/O/N=6 histFitPara = {valA, levelA, wid/8, valB, levelB, wid/8}
FuncFit/NTHR=0/Q DoubleGauss histFitPara histWave /W=histWaveWeight /I=1 /D
if( verbosity > 1 )
ModifyGraph rgb[1] = (0,0,44444)	// fit curve in blue
endif
// step 5: analyze histogram fit	// MAY BE UNNECESSARY
wave fitHistWave = \$("fit_"+histWaveName)
WaveStats/M=1/Q/R=(histFitPara[1],histFitPara[4]) fitHistWave
variable minLevel = V_minLoc, minVal = V_min
if( verbosity > 0 )
printf "threshold level = %g (%g)\r", minLevel, minVal
endif
if( minVal > 0.1 )
print "it is recommended to decimate more since the value at threshold is too high"
endif
// step 6: use FindLevels to evaluate events
FindLevels/Q/DEST=Events sampWave, (levelA+levelB)/2
wave Events
variable numEvents = numpnts(Events)/2
if( verbosity > 0 )
printf "found %g events in total\r", numEvents
endif
Redimension/N=(2, numEvents) Events
SetDataFolder saveDF	// want main output in "current" data folder
Make/D/O/N=(numEvents) EventStart = Events[0][p], LifeTime = Events[1][p] - EventStart
Make/D/O/N=(numEvents) EventLevel = mean(inWave, EventStart, EventStart+LifeTime) - levelA
SetDataFolder pkgDF		// back to our den
// step 7: reject Events of too short duration (lifetime)
variable n
for( n = numEvents - 1; n >= 0; n -= 1 )
DeletePoints n, 1, LifeTime, EventStart, EventLevel
endif
endfor
variable numValidEvents = numpnts(EventStart)
if( verbosity > 0 )
endif
// step 8: display original trace colorized (verbosity > 1)
if( verbosity > 1 )
DoWindow/K DataGraph
Display/N=DataGraph/W=(436.5,42.5,831,251) inWave
string eventWaveName = NameOfWave(inWave)+"_events"
Duplicate/O inWave \$eventWaveName
wave eventWave = \$eventWaveName
variable p1, p2
eventWave = 0
for( n = 0; n < numValidEvents; n += 1 )
p1 = x2pnt(inWave, EventStart[n])
p2 = x2pnt(inWave, EventStart[n] + LifeTime[n])
eventWave[p1,p2] += inWave
endfor
eventWave = eventWave > 0 ? eventWave : NaN
AppendToGraph eventWave
ModifyGraph rgb[1] = (0,0,44444)	// valid events in blue
endif
SetDataFolder saveDF
// KillDataFolder pkgDF
end

//histogram fitting to double gaussian to estimate the step location and stepsize
Function DoubleGauss(w,x) : FitFunc
Wave w
Variable x
return w[0]*exp(-((x-w[1])^2/(2*(w[2]^2))))+w[3]*exp(-((x-w[4])^2/(2*(w[5]^2))))
End```