Merge pull request #1860 from AllenInstitute/feature/1860-expose-psx-…

…deconv-filter-parameters

Expose psx deconv filter parameters

Packages/MIES/MIES_Constants.ipf

@@ -2145,11 +2145,12 @@ Constant SECONDS_PER_DAY = 86400
 Strconstant DB_AXIS_PART_EPOCHS = "_EP"
 /// @}
 
-StrConstant SF_OP_PSX          = "psx"
-StrConstant SF_OP_PSX_KERNEL   = "psxKernel"
-StrConstant SF_OP_PSX_STATS    = "psxStats"
-StrConstant SF_OP_PSX_RISETIME = "psxRiseTime"
-StrConstant SF_OP_PSX_PREP     = "psxPrep"
+StrConstant SF_OP_PSX               = "psx"
+StrConstant SF_OP_PSX_KERNEL        = "psxKernel"
+StrConstant SF_OP_PSX_STATS         = "psxStats"
+StrConstant SF_OP_PSX_RISETIME      = "psxRiseTime"
+StrConstant SF_OP_PSX_PREP          = "psxPrep"
+StrConstant SF_OP_PSX_DECONV_FILTER = "psxDeconvFilter"
 
 /// @name Available PSX states
 /// @anchor PSXStates
@@ -2192,3 +2193,7 @@ StrConstant PSX_STATS_LABELS = "Average;Median;Average Deviation;Standard deviat
 Constant PSX_HORIZ_OFFSET_ONSET = 0
 Constant PSX_HORIZ_OFFSET_PEAK  = 1
 /// @}
+
+Constant PSX_DECONV_FILTER_DEF_LOW   = 0.002
+Constant PSX_DECONV_FILTER_DEF_HIGH  = 0.004
+Constant PSX_DECONV_FILTER_DEF_ORDER = 101

Packages/MIES/MIES_SweepFormula.ipf

@@ -207,7 +207,7 @@ Function/WAVE SF_GetNamedOperations()
 	                  SF_OP_LOG10, SF_OP_APFREQUENCY, SF_OP_CURSORS, SF_OP_SWEEPS, SF_OP_AREA, SF_OP_SETSCALE, SF_OP_BUTTERWORTH, \
 	                  SF_OP_CHANNELS, SF_OP_DATA, SF_OP_LABNOTEBOOK, SF_OP_WAVE, SF_OP_FINDLEVEL, SF_OP_EPOCHS, SF_OP_TP,         \
 	                  SF_OP_STORE, SF_OP_SELECT, SF_OP_POWERSPECTRUM, SF_OP_TPSS, SF_OP_TPBASE, SF_OP_TPINST, SF_OP_TPFIT,        \
-	                  SF_OP_PSX, SF_OP_PSX_KERNEL, SF_OP_PSX_STATS, SF_OP_PSX_RISETIME, SF_OP_PSX_PREP}
+	                  SF_OP_PSX, SF_OP_PSX_KERNEL, SF_OP_PSX_STATS, SF_OP_PSX_RISETIME, SF_OP_PSX_PREP, SF_OP_PSX_DECONV_FILTER}
 
 	return wt
 End
@@ -1099,6 +1099,9 @@ static Function/WAVE SF_FormulaExecutor(string graph, variable jsonID, [string j
 		case SF_OP_PSX_PREP:
 			WAVE out = PSX_OperationPrep(jsonId, jsonPath, graph)
 			break
+		case SF_OP_PSX_DECONV_FILTER:
+			WAVE out = PSX_OperationDeconvFilter(jsonId, jsonPath, graph)
+			break
 		default:
 			SFH_ASSERT(0, "Undefined Operation", jsonId=jsonId)
 	endswitch

Packages/MIES/MIES_SweepFormula_PSX.ipf

@@ -320,9 +320,33 @@ End
 ///
 /// @param sweepData data from a single sweep and channel *without* inserted TP
 /// @param psxKernelFFT FFT'ed kernel from PSX_CreatePSXKernel()
-static Function/WAVE PSX_DeconvoluteSweepData(WAVE sweepData, WAVE/C psxKernelFFT)
+/// @param deconvFilter deconvolution filter settings
+static Function/WAVE PSX_DeconvoluteSweepData(WAVE sweepData, WAVE/C psxKernelFFT, WAVE deconvFilter)
 
-	variable numPoints, fftSize
+	variable numPoints, fftSize, samp, low, high, order, lowFrac, highFrac
+
+	samp  = 1 / (deltax(sweepData) * MILLI_TO_ONE)
+	low   = deconvFilter[%$"Filter Low"]
+	high  = deconvFilter[%$"Filter High"]
+	order = deconvFilter[%$"Filter Order"]
+
+	if(IsNaN(low))
+		lowFrac = PSX_DECONV_FILTER_DEF_LOW
+	else
+		lowFrac = low / samp
+	endif
+
+	if(IsNaN(high))
+		highFrac = PSX_DECONV_FILTER_DEF_HIGH
+	else
+		highFrac = high / samp
+	endif
+
+	if(IsNaN(order))
+		order = PSX_DECONV_FILTER_DEF_ORDER
+	endif
+
+	ASSERT(lowFrac < highFrac, "Expected a low pass filter with lowFrac < highFrac")
 
 	numPoints = DimSize(sweepData, ROWS)
 	fftSize = DimSize(psxKernelFFT, ROWS)
@@ -339,7 +363,7 @@ static Function/WAVE PSX_DeconvoluteSweepData(WAVE sweepData, WAVE/C psxKernelFF
 
 	CopyScales sweepData, Deconv
 
-	FilterFIR/ENDV={3}/LO={0.002, 0.004, 101} Deconv
+	FilterFIR/ENDV={3}/LO={lowFrac, highFrac, order} Deconv
 
 	return Deconv
 end
@@ -406,7 +430,7 @@ end
 /// - deconvolution
 /// - histogram of deconvolution
 /// - gaussian fit of histogram
-static Function [WAVE sweepDataFiltOff, WAVE sweepDataFiltOffDeconv] PSX_Analysis(WAVE sweepData, WAVE psxKernelFFT, variable sweepFilterLow, variable sweepFilterHigh)
+static Function [WAVE sweepDataFiltOff, WAVE sweepDataFiltOffDeconv] PSX_Analysis(WAVE sweepData, WAVE psxKernelFFT, variable sweepFilterLow, variable sweepFilterHigh, WAVE deconvFilter)
 
 	variable offset
 
@@ -419,7 +443,7 @@ static Function [WAVE sweepDataFiltOff, WAVE sweepDataFiltOffDeconv] PSX_Analysi
 		return [$"", $""]
 	endif
 
-	WAVE sweepDataFiltOffDeconv = PSX_DeconvoluteSweepData(sweepDataFiltOff, psxKernelFFT)
+	WAVE sweepDataFiltOffDeconv = PSX_DeconvoluteSweepData(sweepDataFiltOff, psxKernelFFT, deconvFilter)
 
 	return [sweepDataFiltOff, sweepDataFiltOffDeconv]
 end
@@ -696,7 +720,7 @@ static Function/WAVE PSX_CreateOverrideResults(variable numEvents, WAVE/T combos
 	return wv
 End
 
-static Function PSX_OperationSweepGathering(string graph, WAVE/WAVE psxKernelDataset, variable parameterJsonID, variable sweepFilterLow, variable sweepFilterHigh, variable index, WAVE/WAVE output)
+static Function PSX_OperationSweepGathering(string graph, WAVE/WAVE psxKernelDataset, variable parameterJsonID, variable sweepFilterLow, variable sweepFilterHigh, WAVE deconvFilter, variable index, WAVE/WAVE output)
 
 	string key, comboKey, psxParametersAnalyzePeaks, cacheKey
 
@@ -717,7 +741,7 @@ static Function PSX_OperationSweepGathering(string graph, WAVE/WAVE psxKernelDat
 		WAVE sweepDataFiltOff       = psxAnalyzePeaksFromCache[%sweepDataFiltOff]
 		WAVE sweepDataFiltOffDeconv = psxAnalyzePeaksFromCache[%sweepDataFiltOffDeconv]
 	else
-		[WAVE sweepDataFiltOff, WAVE sweepDataFiltOffDeconv] = PSX_Analysis(sweepData, psxKernelFFT, sweepFilterLow, sweepFilterHigh)
+		[WAVE sweepDataFiltOff, WAVE sweepDataFiltOffDeconv] = PSX_Analysis(sweepData, psxKernelFFT, sweepFilterLow, sweepFilterHigh, deconvFilter)
 
 		Make/FREE/WAVE/N=(2) psxAnalyzePeaks
 		SetDimensionLabels(psxAnalyzePeaks, "sweepDataFiltOff;sweepDataFiltOffDeconv", ROWS)
@@ -4099,12 +4123,13 @@ Function/WAVE PSX_Operation(variable jsonId, string jsonPath, string graph)
 	WAVE/WAVE psxKernelDataset = SFH_GetArgumentAsWave(jsonId, jsonPath, graph, SF_OP_PSX, 1,  defOp = "psxKernel()")
 
 	try
-		numberOfSDs   = SFH_GetArgumentAsNumeric(jsonID, jsonPath, graph, SF_OP_PSX, 2, defValue = PSX_NUMBER_OF_SDS_DEFAULT, checkFunc = IsStrictlyPositiveAndFinite)
-		sweepFilterLow     = SFH_GetArgumentAsNumeric(jsonID, jsonPath, graph, SF_OP_PSX, 3, defValue = PSX_DEFAULT_FILTER_LOW)
-		sweepFilterHigh    = SFH_GetArgumentAsNumeric(jsonID, jsonPath, graph, SF_OP_PSX, 4, defValue = PSX_DEFAULT_FILTER_HIGH)
-		maxTauFactor  = SFH_GetArgumentAsNumeric(jsonID, jsonPath, graph, SF_OP_PSX, 5, defValue = PSX_DEFAULT_MAX_TAU_FACTOR, checkFunc = IsStrictlyPositiveAndFinite)
+		numberOfSDs = SFH_GetArgumentAsNumeric(jsonID, jsonPath, graph, SF_OP_PSX, 2, defValue = PSX_NUMBER_OF_SDS_DEFAULT, checkFunc = IsStrictlyPositiveAndFinite)
+		sweepFilterLow = SFH_GetArgumentAsNumeric(jsonID, jsonPath, graph, SF_OP_PSX, 3, defValue = PSX_DEFAULT_FILTER_LOW, checkFunc = IsNullOrPositiveAndFinite)
+		sweepFilterHigh = SFH_GetArgumentAsNumeric(jsonID, jsonPath, graph, SF_OP_PSX, 4, defValue = PSX_DEFAULT_FILTER_HIGH, checkFunc = IsNullOrPositiveAndFinite)
+		maxTauFactor = SFH_GetArgumentAsNumeric(jsonID, jsonPath, graph, SF_OP_PSX, 5, defValue = PSX_DEFAULT_MAX_TAU_FACTOR, checkFunc = IsStrictlyPositiveAndFinite)
 		WAVE riseTime = SFH_GetArgumentAsWave(jsonID, jsonPath, graph, SF_OP_PSX, 6, defOp = "psxRiseTime()", singleResult = 1)
 		ASSERT(IsNumericWave(riseTime), "Invalid return from psxRiseTime")
+		WAVE deconvFilter = SFH_GetArgumentAsWave(jsonID, jsonPath, graph, SF_OP_PSX, 7, defOp = "psxDeconvFilter()", singleResult = 1)
 
 		parameterJsonID = JWN_GetWaveNoteAsJSON(psxKernelDataset)
 		parameterPath = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX
@@ -4118,6 +4143,11 @@ Function/WAVE PSX_Operation(variable jsonId, string jsonPath, string graph)
 		JSON_AddTreeObject(parameterJsonID, parameterPath)
 		JSON_AddVariable(parameterJsonID, parameterPath + "/upperThreshold", riseTime[%$"Upper Threshold"])
 		JSON_AddVariable(parameterJsonID, parameterPath + "/lowerThreshold", riseTime[%$"Lower Threshold"])
+		parameterPath = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX_DECONV_FILTER
+		JSON_AddTreeObject(parameterJsonID, parameterPath)
+		JSON_AddVariable(parameterJsonID, parameterPath + "/filterLow", deconvFilter[%$"Filter Low"])
+		JSON_AddVariable(parameterJsonID, parameterPath + "/filterHigh", deconvFilter[%$"Filter High"])
+		JSON_AddVariable(parameterJsonID, parameterPath + "/filterOrder", deconvFilter[%$"Filter Order"])
 
 		numCombos = DimSize(psxKernelDataset, ROWS) / PSX_KERNEL_OUTPUTWAVES_PER_ENTRY
 		ASSERT(IsInteger(numCombos) && numCombos > 0, "Invalid number of input sets from psxKernel()")
@@ -4138,7 +4168,7 @@ Function/WAVE PSX_Operation(variable jsonId, string jsonPath, string graph)
 		endfor
 
 		for(i = 0; i < numCombos; i += 1)
-			PSX_OperationSweepGathering(graph, psxKernelDataset, parameterJsonID, sweepFilterLow, sweepFilterHigh, i, output)
+			PSX_OperationSweepGathering(graph, psxKernelDataset, parameterJsonID, sweepFilterLow, sweepFilterHigh, deconvFilter, i, output)
 		endfor
 
 		[WAVE hist, WAVE fit, peakThresh, dataUnit] = PSX_CalculatePeakThreshold(output, numCombos, numberOfSDs)
@@ -4253,6 +4283,26 @@ Function/WAVE PSX_OperationRiseTime(variable jsonId, string jsonPath, string gra
 	return SFH_GetOutputForExecutor(output, graph, SF_OP_PSX_RISETIME)
 End
 
+Function/WAVE PSX_OperationDeconvFilter(variable jsonId, string jsonPath, string graph)
+
+	variable low, high, order
+
+	SFH_CheckArgumentCount(jsonId, jsonPath, SF_OP_PSX_DECONV_FILTER, 0, maxArgs = 3)
+
+	low = SFH_GetArgumentAsNumeric(jsonId, jsonPath, graph, SF_OP_PSX_DECONV_FILTER, 0, defValue = NaN, checkFunc = IsNullOrPositiveAndFinite, checkDefault = 0)
+	high = SFH_GetArgumentAsNumeric(jsonId, jsonPath, graph, SF_OP_PSX_DECONV_FILTER, 1, defValue = NaN, checkFunc = IsNullOrPositiveAndFinite, checkDefault = 0)
+	order = SFH_GetArgumentAsNumeric(jsonId, jsonPath, graph, SF_OP_PSX_DECONV_FILTER, 2, defValue = NaN, checkFunc = IsOdd, checkDefault = 0)
+
+	Make/D/FREE params = {low, high, order}
+	SetDimensionLabels(params, "Filter Low;Filter High;Filter Order", ROWS)
+
+	WAVE/WAVE output = SFH_CreateSFRefWave(graph, SF_OP_PSX_DECONV_FILTER, 1)
+
+	output[0] = params
+
+	return SFH_GetOutputForExecutor(output, graph, SF_OP_PSX_DECONV_FILTER)
+End
+
 Function/WAVE PSX_OperationStats(variable jsonId, string jsonPath, string graph)
 
 	string stateAsStr, postProc, id, prop

Packages/doc/SweepFormula.rst

@@ -904,7 +904,7 @@ The `psx` operation allows to classify miniature PSC/PSP's interactively.
 
 .. code-block:: bash
 
-   psx(id, [psxKernel(), numSDs, filterLow, filterHigh, maxTauFactor, psxRiseTime()])
+   psx(id, [psxKernel(), numSDs, filterLow, filterHigh, maxTauFactor, psxRiseTime(), psxDeconvFilter()])
 
 The function accepts one to seven arguments.
 
@@ -931,6 +931,9 @@ maxTauFactor
 psxRiseTime
   results from the `psxRiseTime` operation
 
+psxDeconvFilter
+  results from the `psxDeconvFilter` operation
+
 The plotting is implemented in a custom way. Due to that multiple `psx`
 operations can only be separated by `with` and not `and`.
 
@@ -1015,6 +1018,33 @@ upperThreshold
    psxRiseTime(0.5)
    psxRiseTime(0.5, 0.9)
 
+psxDeconvFilter
+"""""""""""""""
+
+The `psxDeconvFilter` operation is a helper operation for `psx` to manage the deconvolution filter settings.
+
+   psxDeconvFilter([lowFreq, highFreq, order])
+
+The function accepts zero to three arguments.
+
+lowFreq [Hz]
+   defaults to `NaN`
+
+highFreq [Hz]
+   defaults to `NaN`
+
+order
+   defaults to `NaN`
+
+The default values of `NaN` are replaced inside `psx`. For the order this is
+`101`, for the frequencies this is a normalized frequency which depends on the
+sampling interval of the data.
+
+.. code-block:: bash
+
+   psxDeconvFilter(500, 1000)
+   psxDeconvFilter(400, 600, 91)
+
 psxstats
 """"""""
 

Packages/tests/Basic/UTF_SweepFormula_PSX.ipf

@@ -2992,3 +2992,97 @@ static Function TestStoreAndLoad()
 	CheckPSXEventField({psxEvent_0}, {"Fit manual QC call", "Event manual QC call"}, {1}, PSX_UNDET)
 	CheckPSXEventField({psxEvent_1}, {"Fit manual QC call", "Event manual QC call"}, {0}, PSX_UNDET)
 End
+
+static Function [variable filterLow, variable filterHigh, variable filterOrder] TestDevonvFilterContainer(WAVE/WAVE dataWref)
+
+	CHECK_WAVE(dataWref, WAVE_WAVE)
+	CHECK_EQUAL_VAR(DimSize(dataWref, ROWS), 1)
+	WAVE/Z data = dataWref[0]
+	CHECK_WAVE(data, NUMERIC_WAVE)
+	CHECK_EQUAL_VAR(DimSize(data, ROWS), 3)
+
+	filterLow = data[%$"Filter Low"]
+	if(IsNaN(filterLow))
+		PASS()
+	else
+		CHECK(IsNullOrPositiveAndFinite(filterLow))
+	endif
+
+	filterHigh = data[%$"Filter High"]
+	if(IsNaN(filterHigh))
+		PASS()
+	else
+		CHECK(IsNullOrPositiveAndFinite(filterHigh))
+	endif
+
+	filterOrder = data[%$"Filter Order"]
+
+	if(IsNaN(filterOrder))
+		PASS()
+	else
+		CHECK(IsOdd(filterOrder))
+	endif
+
+	return [filterLow, filterHigh, filterOrder]
+End
+
+static Function TestOperationDeconvFilter()
+
+	string win, device, str
+	variable filterLow, filterHigh, filterOrder
+
+	[win, device] = CreateFakeDataBrowserWindow()
+
+	str = "psxDeconvFilter()"
+	WAVE/WAVE dataWref = SF_ExecuteFormula(str, win, useVariables = 0)
+	[filterLow, filterHigh, filterOrder] = TestDevonvFilterContainer(dataWref)
+	CHECK_EQUAL_VAR(filterLow, NaN)
+	CHECK_EQUAL_VAR(filterHigh, NaN)
+	CHECK_EQUAL_VAR(filterOrder, NaN)
+
+	str = "psxDeconvFilter(40)"
+	WAVE/WAVE dataWref = SF_ExecuteFormula(str, win, useVariables = 0)
+	[filterLow, filterHigh, filterOrder] = TestDevonvFilterContainer(dataWref)
+	CHECK_EQUAL_VAR(filterLow, 40)
+	CHECK_EQUAL_VAR(filterHigh, NaN)
+	CHECK_EQUAL_VAR(filterOrder, NaN)
+
+	str = "psxDeconvFilter(40, 50)"
+	WAVE/WAVE dataWref = SF_ExecuteFormula(str, win, useVariables = 0)
+	[filterLow, filterHigh, filterOrder] = TestDevonvFilterContainer(dataWref)
+	CHECK_EQUAL_VAR(filterLow, 40)
+	CHECK_EQUAL_VAR(filterHigh, 50)
+	CHECK_EQUAL_VAR(filterOrder, NaN)
+
+	str = "psxDeconvFilter(40, 50, 11)"
+	WAVE/WAVE dataWref = SF_ExecuteFormula(str, win, useVariables = 0)
+	[filterLow, filterHigh, filterOrder] = TestDevonvFilterContainer(dataWref)
+	CHECK_EQUAL_VAR(filterLow, 40)
+	CHECK_EQUAL_VAR(filterHigh, 50)
+	CHECK_EQUAL_VAR(filterOrder, 11)
+
+	// check parameters
+	try
+		str = "psxDeconvFilter(-1)"
+		WAVE/WAVE dataWref = SF_ExecuteFormula(str, win, useVariables = 0)
+		FAIL()
+	catch
+		CHECK_NO_RTE()
+	endtry
+
+	try
+		str = "psxDeconvFilter(1, -1)"
+		WAVE/WAVE dataWref = SF_ExecuteFormula(str, win, useVariables = 0)
+		FAIL()
+	catch
+		CHECK_NO_RTE()
+	endtry
+
+	try
+		str = "psxDeconvFilter(1, 1, 2)"
+		WAVE/WAVE dataWref = SF_ExecuteFormula(str, win, useVariables = 0)
+		FAIL()
+	catch
+		CHECK_NO_RTE()
+	endtry
+End