diff --git a/Packages/MIES/MIES_Cache.ipf b/Packages/MIES/MIES_Cache.ipf index bda268526f..21c730bdf7 100644 --- a/Packages/MIES/MIES_Cache.ipf +++ b/Packages/MIES/MIES_Cache.ipf @@ -499,7 +499,7 @@ End /// @param psxParameters JSON dump of the psx/psxKernel operation parameters Function/S CA_PSXEventsKey(string comboKey, string psxParameters) - return CA_PSXBaseKey(comboKey, psxParameters) + " Events " + ":Version 2" + return CA_PSXBaseKey(comboKey, psxParameters) + " Events " + ":Version 3" End Function/S CA_PSXOperationKey(string comboKey, string psxParameters) @@ -509,7 +509,7 @@ End Function/S CA_PSXAnalyzePeaks(string comboKey, string psxParameters) - return CA_PSXBaseKey(comboKey, psxParameters) + " Analyze Peaks " + ":Version 2" + return CA_PSXBaseKey(comboKey, psxParameters) + " Analyze Peaks " + ":Version 3" End /// @brief Return the key for the igor info entries @@ -582,6 +582,12 @@ threadsafe Function/S CA_CalculateFetchEpochsKey(WAVE numericalvalues, WAVE text return "Version 1:" + Hash(key + num2istr(crc), HASH_SHA2_256) End + +Function/S CA_GetGoodFFTSizesKeys() + + return "GetGoodFFTSizes Version 1" +End + ///@} /// @brief Make space for one new entry in the cache waves diff --git a/Packages/MIES/MIES_Constants.ipf b/Packages/MIES/MIES_Constants.ipf index a8fb0b502f..6c81526cf2 100644 --- a/Packages/MIES/MIES_Constants.ipf +++ b/Packages/MIES/MIES_Constants.ipf @@ -2358,12 +2358,13 @@ 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_DECONV_FILTER = "psxDeconvFilter" +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_BP_FILTER = "psxDeconvBPFilter" +StrConstant SF_OP_PSX_SWEEP_BP_FILTER = "psxSweepBPFilter" /// @name Available PSX states /// @anchor PSXStates @@ -2392,10 +2393,12 @@ Constant PSX_MARKER_UNDET = 18 /// @name Custom error codes for PSX_FitEventDecay() /// @anchor FitEventDecayCustomErrors ///@{ -Constant PSX_DECAY_FIT_ERROR = -10000 +Constant PSX_DECAY_FIT_ERROR = -10000 +Constant PSX_DECAY_FIT_INVALID_RANGE_ERROR = -10001 + ///@} -StrConstant PSX_STATS_LABELS = "Average;Median;Average Deviation;Standard deviation;Skewness;Kurtosis" +StrConstant PSX_STATS_LABELS = "Average;Median;Average Deviation;Standard deviation;Skewness;Kurtosis;Lower quartile;Upper quartile;Inter-quartile range;Median absolute deviation;Most frequent value" /// @name Horizontal offset modes in all event graph /// @@ -2408,9 +2411,8 @@ Constant PSX_HORIZ_OFFSET_PEAK = 1 Constant PSX_HORIZ_OFFSET_SLEW = 2 ///@} -Constant PSX_DECONV_FILTER_DEF_LOW = 500 -Constant PSX_DECONV_FILTER_DEF_HIGH = 50 -Constant PSX_DECONV_FILTER_DEF_ORDER = 7 +Constant PSX_SWEEP_FILTER_DEF_ORDER = 4 +Constant PSX_DECONV_FILTER_DEF_ORDER = 4 StrConstant PSX_JWN_COMBO_KEYS_NAME = "ComboKeys" diff --git a/Packages/MIES/MIES_MiesUtilities_Algorithm.ipf b/Packages/MIES/MIES_MiesUtilities_Algorithm.ipf index 8ad2263139..93acd628f9 100644 --- a/Packages/MIES/MIES_MiesUtilities_Algorithm.ipf +++ b/Packages/MIES/MIES_MiesUtilities_Algorithm.ipf @@ -641,3 +641,80 @@ threadsafe Function/WAVE FindIndizes(WAVE numericOrTextWave, [variable col, stri return result End + +/// @brief Band‑pass filters a wave while automatically reducing IIR filter +/// order until the output contains no NaNs/Infs and its SEM is not larger than +/// the original (simple ringing detection). +/// +/// @param src – input wave +/// @param fHigh – pass‑band edge frequencies in Hz (Igor’s band‑pass requires fLow > fHigh; the routine swaps them if needed) +/// @param fLow – low part +/// @param maxOrder – starting (maximum) IIR filter order to try (>0) +/// +/// Logic: iteratively lowers the filter order until three conditions are met: +/// 1. FilterIIR executes without error. +/// 2. WaveStats reports V_numNaNs = 0 and V_numInfs = 0. +/// 3. SEM(filtered) ≤ SEM(original). +/// +/// @retval realOrder filter order that finally succeeded (NaN if every order failed) +/// @retval filtered filtered data +Function [variable realOrder, WAVE filtered] BandPassWithRingingDetection(WAVE src, variable fHigh, variable fLow, variable maxOrder) + + variable err, samp, semOrig, offset, i + string msg + + ASSERT(IsInteger(maxOrder) && maxOrder > 0 && isEven(maxOrder), "maxOrder must be a positive and even integer") + + // Igor band-pass expects fLow > fHigh + [fHigh, fLow] = MinMax(fLow, fHigh) + + // Sampling rate (Hz) – assumes X scaling is in milliseconds + samp = 1 / (DeltaX(src) * MILLI_TO_ONE) + + // Pre-compute SEM(original) once + WaveStats/Q src + semOrig = V_sem + offset = v_avg + ASSERT(V_numNaNs == 0 && V_numInfs == 0, "Expected only finite value in input data") + + // Prepare destination wave + Duplicate/FREE src, filtered + + for(i = maxOrder; i > 0; i -= 2) + Multithread filtered = src - offset + + FilterIIR/LO=(fLow / samp)/HI=(fHigh / samp)/DIM=(ROWS)/ORD=(i) filtered; err = GetRTError(1) + if(err) + continue + endif + + WaveStats/Q filtered + if(V_numNaNs > 0 || V_numInfs > 0) + sprintf msg, "V_numNaNs: %g, V_numInfs: %g", V_numNaNs, V_numInfs + DEBUGPRINT(msg) + + // bad numerical output → lower order + continue + endif + + if(V_sem > semOrig) + sprintf msg, "V_sem: %g,semOrig: %g", V_sem, semOrig + DEBUGPRINT(msg) + + // noisier than original → ringing + continue + endif + + Multithread filtered += offset + + sprintf msg, "maxOrder: %g, realOrder: %g", maxOrder, i + DEBUGPRINT(msg) + + return [i, filtered] + endfor + + sprintf msg, "maxOrder: %g, realOrder: NaN", maxOrder + DEBUGPRINT(msg) + + return [NaN, $""] +End diff --git a/Packages/MIES/MIES_SweepFormula.ipf b/Packages/MIES/MIES_SweepFormula.ipf index f4cd3b19eb..29efae66c9 100644 --- a/Packages/MIES/MIES_SweepFormula.ipf +++ b/Packages/MIES/MIES_SweepFormula.ipf @@ -66,7 +66,8 @@ Function/WAVE SF_GetNamedOperations() SF_OP_LOG10, SF_OP_APFREQUENCY, SF_OP_CURSORS, SF_OP_SELECTSWEEPS, SF_OP_AREA, SF_OP_SETSCALE, SF_OP_BUTTERWORTH, \ SF_OP_SELECTCHANNELS, 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_DECONV_FILTER, \ + SF_OP_PSX, SF_OP_PSX_KERNEL, SF_OP_PSX_STATS, SF_OP_PSX_RISETIME, SF_OP_PSX_PREP, SF_OP_PSX_DECONV_BP_FILTER, \ + SF_OP_PSX_SWEEP_BP_FILTER, \ SF_OP_MERGE, SF_OP_FIT, SF_OP_FITLINE, SF_OP_DATASET, SF_OP_SELECTVIS, SF_OP_SELECTCM, SF_OP_SELECTSTIMSET, \ SF_OP_SELECTIVSCCSWEEPQC, SF_OP_SELECTIVSCCSETQC, SF_OP_SELECTRANGE, SF_OP_SELECTEXP, SF_OP_SELECTDEV, \ SF_OP_SELECTEXPANDSCI, SF_OP_SELECTEXPANDRAC, SF_OP_SELECTSETCYCLECOUNT, SF_OP_SELECTSETSWEEPCOUNT, \ diff --git a/Packages/MIES/MIES_SweepFormula_Executor.ipf b/Packages/MIES/MIES_SweepFormula_Executor.ipf index ca0583da6f..4ebc848e8e 100644 --- a/Packages/MIES/MIES_SweepFormula_Executor.ipf +++ b/Packages/MIES/MIES_SweepFormula_Executor.ipf @@ -454,8 +454,11 @@ Function/WAVE SFE_FormulaExecutor(STRUCT SF_ExecutionData &exd, [variable srcLoc case SF_OP_PSX_PREP: WAVE out = PSX_OperationPrep(exdop) break - case SF_OP_PSX_DECONV_FILTER: - WAVE out = PSX_OperationDeconvFilter(exdop) + case SF_OP_PSX_DECONV_BP_FILTER: + WAVE out = PSX_OperationDeconvBPFilter(exdop) + break + case SF_OP_PSX_SWEEP_BP_FILTER: + WAVE out = PSX_OperationSweepBPFilter(exdop) break case SF_OP_MERGE: WAVE out = SFO_OperationMerge(exdop) diff --git a/Packages/MIES/MIES_SweepFormula_PSX.ipf b/Packages/MIES/MIES_SweepFormula_PSX.ipf index b569b76824..01bcd778b5 100644 --- a/Packages/MIES/MIES_SweepFormula_PSX.ipf +++ b/Packages/MIES/MIES_SweepFormula_PSX.ipf @@ -59,14 +59,14 @@ static Constant PSX_KEYBOARD_DIR_LR = 1 static Constant PSX_NUMBER_OF_SDS_DEFAULT = 2.5 -static Constant PSX_TAU_CALC_FACTOR = 2.5 -static Constant PSX_BASELINE_RANGE_FACTOR = 10 +static Constant PSX_TAU_CALC_FACTOR = 1 +static Constant PSX_BASELINE_RANGE_FACTOR = 5 static Constant PSX_FIT_RANGE_FACTOR = 10 static Constant PSX_FIT_RANGE_PERC = 0.9 static Constant PSX_BASELINE_NUM_POINTS_AVERAGE = 5 static Constant PSX_PEAK_RANGE_FACTOR_LEFT = 5 static Constant PSX_PEAK_RANGE_FACTOR_RIGHT = 0.33 -static Constant PSX_PEAK_NUM_HIST_BINS = 20 +static Constant PSX_PEAK_NUM_HIST_BINS = 21 static Constant PSX_NUM_PEAKS_MAX = 2000 @@ -123,7 +123,7 @@ static StrConstant PSX_TUD_COMBO_KEY = "comboKey" static StrConstant PSX_TUD_COMBO_INDEX = "comboIndex" static StrConstant PSX_TUD_BLOCK_INDEX = "blockIndex" -static Constant PSX_GUI_SETTINGS_VERSION = 1 +static Constant PSX_GUI_SETTINGS_VERSION = 4 static StrConstant PSX_GUI_SETTINGS_PSX = "GuiSettingsPSX" @@ -152,6 +152,18 @@ static Constant PSX_CACHE_KEY_ANALYZE_PEAKS = 0x3 static Constant EVENT_INDEX_HORIZONTAL = 0x1 static Constant EVENT_INDEX_VERTICAL = 0x2 +static StrConstant PSX_FREE_AXIS_DECONV = "freeDeconvAxis" +static StrConstant PSX_COLORS_DECONV = "65535,16385,36045" + +static StrConstant PSX_FREE_AXIS_PEAK = "freePeakAxis" +static StrConstant PSX_COLORS_PEAK = "65535,43690,0" + +static StrConstant PSX_FREE_AXIS_BASELINE = "freeBaselineAxis" +static StrConstant PSX_COLORS_BASELINE = "0,0,0" + +static Constant PSX_FILTER_SWEEP = 0x1 +static Constant PSX_FILTER_DECONV = 0x2 + Menu "GraphMarquee" "PSX: Accept Event && Fit", /Q, PSX_MouseEventSelection(PSX_ACCEPT, PSX_STATE_EVENT | PSX_STATE_FIT) "PSX: Reject Event && Fit", /Q, PSX_MouseEventSelection(PSX_REJECT, PSX_STATE_EVENT | PSX_STATE_FIT) @@ -160,6 +172,15 @@ Menu "GraphMarquee" "PSX: Jump to Events", /Q, PSX_JumpToEvents() End +static Function PSX_GetFindPeakBoxSize() + +#ifdef AUTOMATED_TESTING + return 10 +#endif // AUTOMATED_TESTING + + return 10 +End + static Function/S PSX_GetUserDataForWorkingFolder() return PSX_USER_DATA_WORKING_FOLDER @@ -175,7 +196,7 @@ static Function/DF PSX_GetWorkingFolder(string win) return $GetUserData(mainWindow, "", PSX_USER_DATA_WORKING_FOLDER) End -static Function/S PSX_GetSpecialPanel(string win) +Function/S PSX_GetSpecialPanel(string win) return GetMainWindow(win) + "#" + PSX_SPECIAL_EVENT_PANEL End @@ -248,26 +269,34 @@ End /// @brief Filter the sweep data /// /// @param sweepDataOff data from a single sweep and channel *without* inserted TP already offsetted -/// @param high high cutoff [Hz] -/// @param low low cutoff [Hz] -static Function/WAVE PSX_FilterSweepData(WAVE/Z sweepDataOff, variable low, variable high) +/// @param sweepFilter low/high cutoff [Hz] and order +static Function [variable realOrder, WAVE filtered] PSX_FilterSweepData(WAVE/Z sweepDataOff, WAVE sweepFilter) - variable samp, err + variable low, high, maxOrder if(!WaveExists(sweepDataOff)) - return $"" + return [NaN, $""] endif - samp = 1 / (deltax(sweepDataOff) * MILLI_TO_ONE) + low = sweepFilter[%$"Filter Low"] + high = sweepFilter[%$"Filter High"] + maxOrder = sweepFilter[%$"Filter Order"] - ASSERT(low > high, "Expected a band pass filter with low > high") + if(IsNaN(low)) + low = sweepFilter[%$"Filter Low (Default)"] + endif - Duplicate/FREE sweepDataOff, filtered + if(IsNaN(high)) + high = sweepFilter[%$"Filter High (Default)"] + endif - FilterIIR/LO=(low / samp)/HI=(high / samp)/DIM=(ROWS)/ORD=6 filtered; err = GetRTError(1) - SFH_ASSERT(!err, "Error filtering the data, msg: " + GetErrMessage(err)) + if(IsNaN(maxOrder)) + maxOrder = sweepFilter[%$"Filter Order (Default)"] + endif + + [realOrder, WAVE filtered] = BandPassWithRingingDetection(sweepDataOff, high, low, maxOrder) - return filtered + return [realOrder, filtered] End /// @brief Create a histogram of the sweep data @@ -329,45 +358,57 @@ static Function [WAVE sweepDataOff, variable offset] PSX_OffsetSweepData(WAVE sw return [output, offset] End +static Function PSX_AddDefaultFilterFrequencies(WAVE filter, variable riseTau, variable decayTau, variable filterType) + + variable fac + + switch(filterType) + case PSX_FILTER_SWEEP: + fac = 1.5 + break + case PSX_FILTER_DECONV: + fac = 2 + break + default: + FATAL_ERROR("Invalid case") + endswitch + + filter[%$"Filter Low (Default)"] = 1 / (fac * pi * riseTau * MILLI_TO_ONE) + filter[%$"Filter High (Default)"] = 1 / (fac * pi * decayTau * MILLI_TO_ONE) +End + /// @brief Return the deconvoluted sweep data /// /// @param sweepData data from a single sweep and channel *without* inserted TP /// @param psxKernelFFT FFT'ed kernel from PSX_CreatePSXKernel() /// @param deconvFilter deconvolution filter settings -static Function/WAVE PSX_DeconvoluteSweepData(WAVE sweepData, WAVE/C psxKernelFFT, WAVE deconvFilter) +static Function [variable realOrder, WAVE filtered] PSX_DeconvoluteSweepData(WAVE sweepData, WAVE/C psxKernelFFT, WAVE deconvFilter) - variable numPoints, fftSize, samp, low, high, order, lowFrac, highFrac, err + variable numPoints, samp, low, high, maxOrder - samp = 1 / (deltax(sweepData) * MILLI_TO_ONE) - low = deconvFilter[%$"Filter Low"] - high = deconvFilter[%$"Filter High"] - order = deconvFilter[%$"Filter Order"] + low = deconvFilter[%$"Filter Low"] + high = deconvFilter[%$"Filter High"] + maxOrder = deconvFilter[%$"Filter Order"] if(IsNaN(low)) - lowFrac = PSX_DECONV_FILTER_DEF_LOW / samp - else - lowFrac = low / samp + low = deconvFilter[%$"Filter Low (Default)"] endif if(IsNaN(high)) - highFrac = PSX_DECONV_FILTER_DEF_HIGH / samp - else - highFrac = high / samp + high = deconvFilter[%$"Filter High (Default)"] endif - if(IsNaN(order)) - order = PSX_DECONV_FILTER_DEF_ORDER + if(IsNaN(maxOrder)) + maxOrder = deconvFilter[%$"Filter Order (Default)"] endif - ASSERT(lowFrac > highFrac, "Expected a band pass filter with low > high") - numPoints = DimSize(sweepData, ROWS) - fftSize = DimSize(psxKernelFFT, ROWS) // no window function on purpose WAVE/C outputFFT = DoFFT(sweepData, padSize = numPoints) - Multithread outputFFT[] = outputFFT[p] / psxKernelFFT[p] + ASSERT(DimSize(outputFFT, ROWS) == DimSize(psxKernelFFT, ROWS), "Unmatched wave sizes") + Multithread outputFFT[] = outputFFT[p] / (psxKernelFFT[p] + 1e-5) IFFT/DEST=Deconv/FREE outputFFT @@ -375,14 +416,9 @@ static Function/WAVE PSX_DeconvoluteSweepData(WAVE sweepData, WAVE/C psxKernelFF ASSERT(V_Value == -1, "Can not handle NaN in the deconvoluted wave") CopyScales sweepData, Deconv - FilterIIR/LO=(lowFrac)/HI=(highFrac)/ORD=(order) Deconv; err = GetRTError(0) + [realOrder, filtered] = BandPassWithRingingDetection(Deconv, high, low, maxOrder) - if(err) - printf "Error applying deconvolution filter: %s\r", GetRTErrMessage() - ClearRTError() - endif - - return Deconv + return [realOrder, filtered] End /// @brief Creates a histogram of the deconvoluted sweep data @@ -430,21 +466,21 @@ End /// - deconvolution /// - histogram of deconvolution /// - gaussian fit of histogram -static Function [WAVE sweepDataOffFilt, WAVE sweepDataOffFiltDeconv] PSX_Analysis(WAVE sweepData, WAVE psxKernelFFT, variable sweepFilterLow, variable sweepFilterHigh, WAVE deconvFilter) +static Function [variable realOrderSweep, variable realOrderDeconv, WAVE sweepDataOffFilt, WAVE sweepDataOffFiltDeconv] PSX_Analysis(WAVE sweepData, WAVE psxKernelFFT, WAVE sweepFilter, WAVE deconvFilter) variable offset [WAVE sweepDataOff, offset] = PSX_OffsetSweepData(sweepData) - WAVE sweepDataOffFilt = PSX_FilterSweepData(sweepDataOff, sweepFilterLow, sweepFilterHigh) + [realOrderSweep, WAVE sweepDataOffFilt] = PSX_FilterSweepData(sweepDataOff, sweepFilter) if(!WaveExists(sweepDataOffFilt)) - return [$"", $""] + return [NaN, NaN, $"", $""] endif - WAVE sweepDataOffFiltDeconv = PSX_DeconvoluteSweepData(sweepDataOffFilt, psxKernelFFT, deconvFilter) + [realOrderDeconv, WAVE sweepDataOffFiltDeconv] = PSX_DeconvoluteSweepData(sweepDataOffFilt, psxKernelFFT, deconvFilter) - return [sweepDataOffFilt, sweepDataOffFiltDeconv] + return [realOrderSweep, realOrderDeconv, sweepDataOffFilt, sweepDataOffFiltDeconv] End /// Searches for peaks in sweepData @@ -459,7 +495,7 @@ End /// @retval peakY y-coordinates of peaks static Function [WAVE/D peakX, WAVE/D peakY] PSX_FindPeaks(WAVE sweepDataOffFiltDeconv, variable threshold, [variable numPeaksMax, variable start, variable stop]) - variable i + variable i, boxSize if(ParamIsDefault(numPeaksMax)) numPeaksMax = PSX_NUM_PEAKS_MAX @@ -473,10 +509,12 @@ static Function [WAVE/D peakX, WAVE/D peakY] PSX_FindPeaks(WAVE sweepDataOffFilt stop = rightx(sweepDataOffFiltDeconv) endif + boxSize = PSX_GetFindPeakBoxSize() + Make/FREE/D/N=(numPeaksMax) peakX, peakY for(i = 0; i < numPeaksMax; i += 1) - FindPeak/B=10/M=(threshold)/Q/R=(start, stop) sweepDataOffFiltDeconv + FindPeak/B=(boxSize)/M=(threshold)/Q/R=(start, stop) sweepDataOffFiltDeconv if(V_Flag != 0) break @@ -505,6 +543,7 @@ static Function [WAVE/D peakX, WAVE/D peakY] PSX_FilterEventsKernelAmpSign(WAVE/ variable peak, peak_t, baseline, baseline_t, amplitude variable overrideSignQC = NaN string comboKey + variable peak_t_prev = -Inf if(!WaveExists(peakXUnfiltered) || !WaveExists(peakYUnfiltered)) return [$"", $""] @@ -516,7 +555,7 @@ static Function [WAVE/D peakX, WAVE/D peakY] PSX_FilterEventsKernelAmpSign(WAVE/ for(i = 0; i < numCrossings; i += 1) - [peak, peak_t, baseline, baseline_t, amplitude] = PSX_CalculateEventProperties(peakXUnfiltered, peakYUnfiltered, sweepDataOffFilt, kernelAmp, kernelRiseTau, kernelDecayTau, i) + [peak, peak_t, baseline, baseline_t, amplitude] = PSX_CalculateEventProperties(peakXUnfiltered, peakYUnfiltered, sweepDataOffFilt, peak_t_prev, kernelAmp, kernelRiseTau, kernelDecayTau, i) #ifdef AUTOMATED_TESTING WAVE/Z overrideResults = GetOverrideResults() @@ -536,10 +575,13 @@ static Function [WAVE/D peakX, WAVE/D peakY] PSX_FilterEventsKernelAmpSign(WAVE/ continue endif - peakX[idx] = peakXUnfiltered[i] - peakY[idx] = peakYUnfiltered[i] + if(isFinite(peakXUnfiltered[i]) && isFinite(peakYUnfiltered[i])) + peakX[idx] = peakXUnfiltered[i] + peakY[idx] = peakYUnfiltered[i] - idx += 1 + idx += 1 + peak_t_prev = peak_t + endif endfor if(idx == 0) @@ -551,23 +593,13 @@ static Function [WAVE/D peakX, WAVE/D peakY] PSX_FilterEventsKernelAmpSign(WAVE/ return [peakX, peakY] End -static Function [variable peak_t, variable peak] PSX_CalculateEventPeak(WAVE peakX, WAVE peakY, WAVE sweepDataOffFilt, variable kernelAmp, variable kernelRiseTau, variable kernelDecayTau, variable index) +static Function [variable peak_t, variable peak] PSX_CalculateEventPeak(WAVE peakX, WAVE peakY, WAVE sweepDataOffFilt, variable baseline_t, variable kernelAmp, variable kernelRiseTau, variable kernelDecayTau, variable index, variable peak_t_prev, variable prevDeconvPeak_t) - variable numCrossings, deconvPeak_t, prevDeconvPeak_t, nextDeconvPeak_t - variable peak_start_search, peak_end_search - - numCrossings = DimSize(peakX, ROWS) - - deconvPeak_t = peakX[index] - - // lower bound - if(index > 0) - prevDeconvPeak_t = peakX[index - 1] - else - prevDeconvPeak_t = -Inf - endif + variable numCrossings, nextDeconvPeak_t + variable peak_start_search, peak_end_search, threshold, range - peak_start_search = max(deconvPeak_t - PSX_PEAK_RANGE_FACTOR_LEFT * kernelRiseTau, prevDeconvPeak_t) + numCrossings = DimSize(peakX, ROWS) + peak_end_search = max(peakX[index] - PSX_PEAK_RANGE_FACTOR_LEFT * kernelRiseTau, baseline_t, prevDeconvPeak_t, peak_t_prev) // upper bound if(index < (numCrossings - 1)) @@ -576,48 +608,88 @@ static Function [variable peak_t, variable peak] PSX_CalculateEventPeak(WAVE pea nextDeconvPeak_t = Inf endif - peak_end_search = min(deconvPeak_t + PSX_PEAK_RANGE_FACTOR_RIGHT * kernelDecayTau, nextDeconvPeak_t) + peak_start_search = min(peakX[index] + PSX_PEAK_RANGE_FACTOR_RIGHT * kernelDecayTau, nextDeconvPeak_t) WaveStats/M=1/Q/R=(peak_start_search, peak_end_search) sweepDataOffFilt - + variable boxSize = 10 + variable ampFactor = 0.25 if(kernelAmp > 0) - peak = V_max - peak_t = V_maxloc + threshold = sweepDataOffFilt(baseline_t) + ampFactor * (V_max - sweepDataOffFilt(baseline_t)) + FindPeak/B=(boxSize)/M=(threshold)/Q/R=(peak_start_search, peak_end_search) sweepDataOffFilt + + if(!V_flag) + peak_t = V_PeakLoc + peak = V_PeakVal + else + range = 2 * DimDelta(sweepDataOffFilt, ROWS) + WaveStats/M=1/Q/R=(peakX[index] - range, peakX[index] + range) sweepDataOffFilt + peak_t = peakX[index] + peak = V_avg + endif elseif(kernelAmp < 0) - peak = V_min - peak_t = V_minloc + threshold = sweepDataOffFilt(baseline_t) + ampFactor * (V_min - sweepDataOffFilt(baseline_t)) + FindPeak/B=(boxSize)/N/M=(threshold)/Q/R=(peak_start_search, peak_end_search) sweepDataOffFilt + + if(!V_flag) + peak_t = V_PeakLoc + peak = V_PeakVal + else + range = 2 * DimDelta(sweepDataOffFilt, ROWS) + WaveStats/M=1/Q/R=(peakX[index] - range, peakX[index] + range) sweepDataOffFilt + peak_t = peakX[index] + peak = V_avg + endif else FATAL_ERROR("Can't handle kernelAmp of zero") endif + ASSERT(IsFinite(peak_t), "peak_t is not finite") + ASSERT(IsFinite(peak), "peak is not finite") + return [peak_t, peak] End -static Function [variable baseline_t, variable baseline] PSX_CalculateEventBaseline(WAVE sweepDataOffFilt, variable peak_t, variable kernelAmp, variable kernelRiseTau) +static Function [variable baseline_t, variable baseline] PSX_CalculateEventBaseline(WAVE sweepDataOffFilt, variable peak_t_prev, variable deconvPeak_t, variable kernelAmp, variable kernelRiseTau, variable prevDeconvPeak_t) + + variable range, start + string str - variable range + start = max(deconvPeak_t - PSX_BASELINE_RANGE_FACTOR * kernelRiseTau, min(peak_t_prev, prevDeconvPeak_t)) - WaveStats/M=1/Q/R=(peak_t - PSX_BASELINE_RANGE_FACTOR * kernelRiseTau, peak_t) sweepDataOffFilt + // the end of the search window for the baseline time is the deconvolved peak time of the current event + WaveStats/M=1/Q/R=(start, deconvPeak_t) sweepDataOffFilt if(kernelAmp > 0) baseline_t = V_minloc elseif(kernelAmp < 0) baseline_t = V_maxloc else - FATAL_ERROR("Can't handle kernelAmp of zero") + sprintf str, "Can't handle kernelAmp of: %g\r", kernelAmp + FATAL_ERROR(str) endif range = PSX_BASELINE_NUM_POINTS_AVERAGE * DimDelta(sweepDataOffFilt, ROWS) WaveStats/M=1/Q/R=(baseline_t - range, baseline_t + range) sweepDataOffFilt baseline = V_avg + ASSERT(isFinite(baseline), "basline is not finite") + return [baseline_t, baseline] End -static Function [variable peak, variable peak_t, variable baseline, variable baseline_t, variable amplitude] PSX_CalculateEventProperties(WAVE peakX, WAVE peakY, WAVE sweepDataOffFilt, variable kernelAmp, variable kernelRiseTau, variable kernelDecayTau, variable index) +static Function [variable peak, variable peak_t, variable baseline, variable baseline_t, variable amplitude] PSX_CalculateEventProperties(WAVE peakX, WAVE peakY, WAVE sweepDataOffFilt, variable peak_t_prev, variable kernelAmp, variable kernelRiseTau, variable kernelDecayTau, variable index) - [peak_t, peak] = PSX_CalculateEventPeak(peakX, peakY, sweepDataOffFilt, kernelAmp, kernelRiseTau, kernelDecayTau, index) - [baseline_t, baseline] = PSX_CalculateEventBaseline(sweepDataOffFilt, peak_t, kernelAmp, kernelRiseTau) + variable deconvPeak_t, prevdeconvPeak_t + + deconvPeak_t = peakX[index] + if(index == 0) + prevdeconvPeak_t = 0 + else + prevdeconvPeak_t = peakX[index - 1] + endif + + [baseline_t, baseline] = PSX_CalculateEventBaseline(sweepDataOffFilt, peak_t_prev, deconvPeak_t, kernelAmp, kernelRiseTau, prevdeconvPeak_t) + [peak_t, peak] = PSX_CalculateEventPeak(peakX, peakY, sweepDataOffFilt, baseline_t, kernelAmp, kernelRiseTau, kernelDecayTau, index, peak_t_prev, prevdeconvPeak_t) amplitude = peak - baseline @@ -628,6 +700,7 @@ End static Function [WAVE/D peakX, WAVE/D peakY] PSX_AnalyzePeaks(WAVE sweepDataOffFiltDeconv, WAVE sweepDataOffFilt, WAVE sweepData, WAVE/Z peakXUnfiltered, WAVE/Z peakYUnfiltered, variable maxTauFactor, variable kernelAmp, variable kernelRiseTau, variable kernelDecayTau, WAVE riseTimeParams, WAVE psxEvent, WAVE eventFit) variable i, numCrossings, deconvPeak, deconvPeak_t, peak, peak_t, baseline, baseline_t, amplitude, iei + variable peak_t_prev = -Inf // sets the previous peak time to the start of the range // we need to first throw away events with invalid amplitude so that // we can then calculate the distance to the neighbour in peakX[i + 1] below @@ -649,8 +722,8 @@ static Function [WAVE/D peakX, WAVE/D peakY] PSX_AnalyzePeaks(WAVE sweepDataOffF deconvPeak_t = peakX[i] deconvPeak = peakY[i] - [peak, peak_t, baseline, baseline_t, amplitude] = PSX_CalculateEventProperties(peakX, peakY, sweepDataOffFilt, \ - kernelAmp, kernelRiseTau, kernelDecayTau, i) + [peak, peak_t, baseline, baseline_t, amplitude] = PSX_CalculateEventProperties(peakX, peakY, sweepDataOffFilt, \ + peak_t_prev, kernelAmp, kernelRiseTau, kernelDecayTau, i) if(i == 0) iei = NaN @@ -667,6 +740,8 @@ static Function [WAVE/D peakX, WAVE/D peakY] PSX_AnalyzePeaks(WAVE sweepDataOffF psxEvent[i][%baseline_t] = baseline_t psxEvent[i][%amplitude] = amplitude psxEvent[i][%iei] = iei + + peak_t_prev = peak_t endfor // safe defaults @@ -682,12 +757,17 @@ static Function [WAVE/D peakX, WAVE/D peakY] PSX_AnalyzePeaks(WAVE sweepDataOffF riseTimeParams[%$"Differentiate Threshold"], \ p) - Multithread psxEvent[][%$"Rise Time"] = PSX_CalculateRiseTime(sweepDataOffFilt, psxEvent, kernelAmp, \ - riseTimeParams[%$"Lower Threshold"], \ - riseTimeParams[%$"Upper Threshold"], \ - p) + Multithread psxEvent[][%$"Onset"] = GetInterpolatedYValue(sweepDataOffFilt, psxEvent[p][%$"Onset Time"]) - psxEvent[][%tau] = PSX_FitEventDecay(sweepDataOffFilt, psxEvent, maxTauFactor, eventFit, p) + Multithread psxEvent[][%$"Rise Time"] = PSX_CalculateRiseTimeWrapper(sweepDataOffFilt, psxEvent, kernelAmp, \ + riseTimeParams[%$"Lower Threshold"], \ + riseTimeParams[%$"Upper Threshold"], \ + p) + + Multithread psxEvent[][%$"Rise"] = GetInterpolatedYValue(sweepDataOffFilt, psxEvent[p][%$"Rise Time"]) + + Make/FREE/N=(DimSize(psxEvent, ROWS)) indexHelper + indexHelper[] = PSX_FitEventDecay(sweepDataOffFilt, psxEvent, maxTauFactor, eventFit, p) return [peakX, peakY] End @@ -712,20 +792,26 @@ End // @brief Returns a good tau which does capture a lot of the tau events static Function PSX_GetGoodTauImpl(WAVE psxEvent) - variable numEvents, err, xVal, idx + variable numEvents, err, xVal, idx, tau idx = FindDimLabel(psxEvent, COLS, "tau") Duplicate/FREE/RMD=[][idx] psxEvent, tauWithNaN - WAVE/Z tau = ZapNaNs(tauWithNaN) + WAVE/Z taus = ZapNaNs(tauWithNaN) - if(!WaveExists(tau)) + if(!WaveExists(taus)) return PSX_DEFAULT_X_START_OFFSET endif - WaveStats/Q tau + WaveStats/Q taus - return V_avg + PSX_TAU_CALC_FACTOR * V_sdev + tau = V_avg + PSX_TAU_CALC_FACTOR * V_sdev + + if(IsFinite(tau)) + return tau + endif + + return PSX_DEFAULT_X_START_OFFSET End /// @brief Return the x-axis range useful for displaying and extracting a single event @@ -754,6 +840,8 @@ static Function [variable first, variable last] PSX_GetSingleEventRange(WAVE psx last = min(first + offset, psxEvent[index + 1][%baseline_t]) endif + ASSERT(first <= last, "range must have first < last") + return [first, last] End @@ -762,25 +850,25 @@ End /// x-zero is taken from sweepData static Function [variable start, variable stop] PSX_GetEventFitRange(WAVE sweepDataOffFilt, WAVE psxEvent, variable eventIndex) - variable calcLength, maxLength + variable calcLength, maxLength, decayTau - start = psxEvent[eventIndex][%deconvPeak_t] + start = psxEvent[eventIndex][%peak_t] - maxLength = PSX_FIT_RANGE_FACTOR * JWN_GetNumberFromWaveNote(psxEvent, SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/psxKernel/decayTau") + decayTau = JWN_GetNumberFromWaveNote(psxEvent, SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/psxKernel/decayTau") + maxLength = PSX_FIT_RANGE_FACTOR * decayTau + ASSERT(isFinite(maxLength), "Failed to retrieve finite decay tau") if(eventIndex == (DimSize(psxEvent, ROWS) - 1)) - calcLength = maxLength + stop = start + maxLength else - calcLength = min((psxEvent[eventIndex + 1][%deconvPeak_t] - start) * PSX_FIT_RANGE_PERC, maxLength) - endif - - if(calcLength == 0) - calcLength = maxLength + stop = min(start + 5 * decayTau, psxEvent[eventIndex + 1][%baseline_t]) // min(psxEvent[eventIndex + 1][%peak_t] * PSX_FIT_RANGE_PERC, psxEvent[eventIndex + 1][%baseline_t]) endif - stop = min(start + calcLength, IndexToScale(sweepDataOffFilt, DimSize(sweepDataOffFilt, ROWS), ROWS)) + stop = min(stop, IndexToScale(sweepDataOffFilt, DimSize(sweepDataOffFilt, ROWS), ROWS)) - ASSERT(start < stop, "Invalid fit range calculation") + if(start > stop) + return [NaN, NaN] + endif return [start, stop] End @@ -790,32 +878,43 @@ End /// /// \rst /// -/// exp_XOffset: :math:`y = K0 + K1 \cdot exp(-(x - x0)/K2)` +/// dblexp_XOffset: :math:`y = K0 + K1 \cdot exp(-(x - x0) / K2) + K3 * exp(-(x - x0) / K4)` /// /// \endrst static Function PSX_FitEventDecay(WAVE sweepDataOffFilt, WAVE psxEvent, variable maxTauFactor, WAVE/WAVE eventFit, variable eventIndex) - variable startTime, endTime, err, decayTau, fitRange, overrideTau + variable startTime, endTime, err, weightedTau, slowTau, fastTau, fitRange, overrideTau string comboKey [startTime, endTime] = PSX_GetEventFitRange(sweepDataOffFilt, psxEvent, eventIndex) + if(IsNaN(startTime) && IsNaN(endTime)) + psxEvent[eventIndex][%$"Fit manual QC call"] = PSX_REJECT + psxEvent[eventIndex][%$"Fit result"] = PSX_DECAY_FIT_INVALID_RANGE_ERROR + psxEvent[eventIndex][%weightedTau] = NaN + psxEvent[eventIndex][%slowTau] = NaN + psxEvent[eventIndex][%fastTau] = NaN + + return NaN + endif + DFREF currDFR = GetDataFolderDFR() SetDataFolder NewFreeDataFolder() // require a converging exponential Make/FREE/T constraints = {"K2 > 0"} - Make/FREE/D/N=3 coefWave + Make/FREE/D/N=5 coefWave AssertOnAndClearRTError() - CurveFit/Q/N=1/NTHR=1/M=0/W=2 exp_XOffset, kwCWave=coefWave, sweepDataOffFilt(startTime, endTime)/D/C=constraints; err = GetRTError(1) + CurveFit/Q/N=1/NTHR=1/M=0/W=2 dblexp_XOffset, kwCWave=coefWave, sweepDataOffFilt(startTime, endTime)/D/C=constraints; err = GetRTError(1) WAVE fit = MakeWaveFree($"fit__free_") SetDataFolder currDFR - - decayTau = coefWave[2] + // weighted tau computed from double exponential fit + [fastTau, slowTau] = MinMax(coefWave[2], coefWave[4]) + weightedTau = ((coefWave[1] * slowTau + coefWave[3] * fastTau) / (coefWave[1] + coefWave[3])) #ifdef AUTOMATED_TESTING WAVE/Z overrideResults = GetOverrideResults() @@ -829,7 +928,7 @@ static Function PSX_FitEventDecay(WAVE sweepDataOffFilt, WAVE psxEvent, variable overrideTau = overrideResults[eventIndex][%$comboKey][%Tau] if(!IsNaN(overrideTau)) - decayTau = overrideTau + weightedTau = overrideTau endif endif #endif // AUTOMATED_TESTING @@ -837,14 +936,22 @@ static Function PSX_FitEventDecay(WAVE sweepDataOffFilt, WAVE psxEvent, variable if(err) psxEvent[eventIndex][%$"Fit manual QC call"] = PSX_REJECT psxEvent[eventIndex][%$"Fit result"] = -err + psxEvent[eventIndex][%weightedTau] = NaN + psxEvent[eventIndex][%slowTau] = NaN + psxEvent[eventIndex][%fastTau] = NaN + return NaN endif fitRange = endTime - startTime - if(IsFinite(decayTau) && decayTau > (maxTauFactor * fitRange)) + if((IsFinite(weightedTau) && weightedTau > (maxTauFactor * fitRange)) || fastTau <= 0 || slowTau <= 0 || weightedTau <= 0) psxEvent[eventIndex][%$"Fit manual QC call"] = PSX_REJECT psxEvent[eventIndex][%$"Fit result"] = PSX_DECAY_FIT_ERROR + psxEvent[eventIndex][%weightedTau] = NaN + psxEvent[eventIndex][%slowTau] = NaN + psxEvent[eventIndex][%fastTau] = NaN + return NaN endif @@ -854,8 +961,9 @@ static Function PSX_FitEventDecay(WAVE sweepDataOffFilt, WAVE psxEvent, variable eventFit[eventIndex] = fit psxEvent[eventIndex][%$"Fit result"] = 1 psxEvent[eventIndex][%$"Fit manual QC call"] = PSX_UNDET - - return decayTau + psxEvent[eventIndex][%weightedTau] = weightedTau + psxEvent[eventIndex][%slowTau] = slowTau + psxEvent[eventIndex][%fastTau] = fastTau End /// @brief Create the override results 2D wave @@ -888,14 +996,15 @@ static Function/WAVE PSX_CreateOverrideResults(variable numEvents, WAVE/T combos End /// @return 0 on success, 1 otherwise -static Function PSX_OperationSweepGathering(string graph, WAVE/WAVE psxKernelDataset, variable parameterJsonID, variable sweepFilterLow, variable sweepFilterHigh, WAVE deconvFilter, variable index, WAVE/WAVE output) +static Function PSX_OperationSweepGathering(string graph, WAVE/WAVE psxKernelDataset, variable parameterJsonID, WAVE sweepFilter, WAVE deconvFilter, variable readIndex, variable writeIndex, WAVE/WAVE output) string key, comboKey, psxParametersAnalyzePeaks, cacheKey + variable realOrderSweep, realOrderDeconv - key = PSX_GenerateKey("psxKernelFFT", index) + key = PSX_GenerateKey("psxKernelFFT", readIndex) WAVE psxKernelFFT = psxKernelDataset[%$key] - key = PSX_GenerateKey("sweepData", index) + key = PSX_GenerateKey("sweepData", readIndex) WAVE sweepData = psxKernelDataset[%$key] [WAVE selectData, WAVE range] = SFH_ParseToSelectDataWaveAndRange(sweepData) @@ -914,16 +1023,24 @@ static Function PSX_OperationSweepGathering(string graph, WAVE/WAVE psxKernelDat WAVE sweepDataOffFilt = psxAnalyzePeaksFromCache[%sweepDataOffFilt] WAVE sweepDataOffFiltDeconv = psxAnalyzePeaksFromCache[%sweepDataOffFiltDeconv] + realOrderSweep = WaveRef(psxAnalyzePeaksFromCache[%realOrderSweep])[0] + realOrderDeconv = WaveRef(psxAnalyzePeaksFromCache[%realOrderDeconv])[0] else - [WAVE sweepDataOffFilt, WAVE sweepDataOffFiltDeconv] = PSX_Analysis(sweepData, psxKernelFFT, sweepFilterLow, sweepFilterHigh, deconvFilter) + [realOrderSweep, realOrderDeconv, WAVE sweepDataOffFilt, WAVE sweepDataOffFiltDeconv] = PSX_Analysis(sweepData, psxKernelFFT, sweepFilter, deconvFilter) if(!WaveExists(sweepDataOffFilt) || !WaveExists(sweepDataOffFiltDeconv)) Make/FREE/WAVE/N=(0) psxAnalyzePeaks else - Make/FREE/WAVE/N=(2) psxAnalyzePeaks - SetDimensionLabels(psxAnalyzePeaks, "sweepDataOffFilt;sweepDataOffFiltDeconv", ROWS) + Make/FREE/WAVE/N=(4) psxAnalyzePeaks + SetDimensionLabels(psxAnalyzePeaks, "sweepDataOffFilt;sweepDataOffFiltDeconv;realOrderSweep;realOrderDeconv", ROWS) psxAnalyzePeaks[%sweepDataOffFilt] = sweepDataOffFilt psxAnalyzePeaks[%sweepDataOffFiltDeconv] = sweepDataOffFiltDeconv + + Make/FREE pkg = {realOrderSweep} + psxAnalyzePeaks[%realOrderSweep] = pkg + + Make/FREE pkg = {realOrderDeconv} + psxAnalyzePeaks[%realOrderDeconv] = pkg endif CA_StoreEntryIntoCache(cacheKey, psxAnalyzePeaks) @@ -933,13 +1050,17 @@ static Function PSX_OperationSweepGathering(string graph, WAVE/WAVE psxKernelDat endif endif - key = PSX_GenerateKey("sweepData", index) + key = PSX_GenerateKey("sweepData", writeIndex) output[%$key] = sweepData - key = PSX_GenerateKey("sweepDataOffFilt", index) + JWN_SetNumberInWaveNote(sweepDataOffFilt, "maxOrder", sweepFilter[%$"Filter Order"]) + JWN_SetNumberInWaveNote(sweepDataOffFilt, "realOrder", realOrderSweep) + key = PSX_GenerateKey("sweepDataOffFilt", writeIndex) output[%$key] = sweepDataOffFilt - key = PSX_GenerateKey("sweepDataOffFiltDeconv", index) + JWN_SetNumberInWaveNote(sweepDataOffFiltDeconv, "maxOrder", deconvFilter[%$"Filter Order"]) + JWN_SetNumberInWaveNote(sweepDataOffFiltDeconv, "realOrder", realOrderDeconv) + key = PSX_GenerateKey("sweepDataOffFiltDeconv", writeIndex) output[%$key] = sweepDataOffFiltDeconv return 0 @@ -1145,8 +1266,14 @@ static Function [WAVE/D results, WAVE eventIndex, WAVE marker, WAVE/T comboKeys] case "xinterval": propLabel = "iei" break - case "tau": - propLabel = "tau" + case "slowtau": + propLabel = "slowTau" + break + case "fasttau": + propLabel = "fastTau" + break + case "weightedtau": + propLabel = "weightedTau" break case "estate": propLabel = "Event manual QC call" @@ -1163,11 +1290,17 @@ static Function [WAVE/D results, WAVE eventIndex, WAVE marker, WAVE/T comboKeys] case "slewratetime": propLabel = "Slew Rate Time" break + case "onsettime": + propLabel = "Onset Time" + break + case "onset": + propLabel = "Onset" + break case "risetime": propLabel = "Rise Time" break - case "onsettime": - propLabel = "Onset Time" + case "rise": + propLabel = "Rise" break default: FATAL_ERROR("Impossible prop: " + prop) @@ -1191,7 +1324,9 @@ static Function [WAVE/D results, WAVE eventIndex, WAVE marker, WAVE/T comboKeys] stateType = "Event manual QC call" break case "Fit result": // fallthrough - case "tau": // fallthrough + case "slowtau": // fallthrough + case "fasttau": // fallthrough + case "weightedtau": // fallthrough case "Fit manual QC call": stateType = "Fit manual QC call" break @@ -1485,8 +1620,14 @@ static Function/WAVE PSX_OperationStatsImpl(string graph, string id, WAVE/WAVE s case "xinterval": propLabelAxis = "Event interval" + " (" + JWN_GetStringFromWaveNote(allEvents[0], PSX_X_DATA_UNIT) + ")" break - case "tau": - propLabelAxis = "Decay tau" + " (" + JWN_GetStringFromWaveNote(allEvents[0], PSX_X_DATA_UNIT) + ")" + case "slowtau": + propLabelAxis = "Slow tau" + " (" + JWN_GetStringFromWaveNote(allEvents[0], PSX_X_DATA_UNIT) + ")" + break + case "fasttau": + propLabelAxis = "Fast tau" + " (" + JWN_GetStringFromWaveNote(allEvents[0], PSX_X_DATA_UNIT) + ")" + break + case "weightedtau": + propLabelAxis = "Weighted tau" + " (" + JWN_GetStringFromWaveNote(allEvents[0], PSX_X_DATA_UNIT) + ")" break case "estate": propLabelAxis = "Event manual QC" + " (enum)" @@ -1503,13 +1644,18 @@ static Function/WAVE PSX_OperationStatsImpl(string graph, string id, WAVE/WAVE s case "slewratetime": propLabelAxis = "Slew Rate time" + " (" + JWN_GetStringFromWaveNote(allEvents[0], PSX_X_DATA_UNIT) + ")" break + case "onsettime": + propLabelAxis = "Onset time" + " (" + JWN_GetStringFromWaveNote(allEvents[0], PSX_X_DATA_UNIT) + ")" + break + case "onset": + propLabelAxis = "Onset" + " (" + JWN_GetStringFromWaveNote(allEvents[0], PSX_Y_DATA_UNIT) + ")" + break case "risetime": propLabelAxis = "Rise time" + " (" + JWN_GetStringFromWaveNote(allEvents[0], PSX_X_DATA_UNIT) + ")" break - case "onsettime": - propLabelAxis = "Onset time" + " (" + JWN_GetStringFromWaveNote(allEvents[0], PSX_X_DATA_UNIT) + ")" + case "rise": + propLabelAxis = "Rise" + " (" + JWN_GetStringFromWaveNote(allEvents[0], PSX_Y_DATA_UNIT) + ")" break - default: FATAL_ERROR("Impossible prop: " + prop) endswitch @@ -1543,13 +1689,19 @@ static Function/WAVE PSX_OperationStatsImpl(string graph, string id, WAVE/WAVE s WaveStats/Q/M=2 resultsRawClean - Make/FREE/D results = {V_avg, NaN, V_adev, V_sdev, V_skew, V_kurt} + Make/FREE/D results = {V_avg, NaN, V_adev, V_sdev, V_skew, V_kurt, NaN, NaN, NaN, NaN, NaN} + SetDimensionLabels(results, PSX_STATS_LABELS, ROWS) StatsQuantiles/Q/Z resultsRawClean MakeWaveFree($"W_StatsQuantiles") if(!V_Flag) - results[1] = V_Median + results[%$"Median"] = V_Median + results[%$"Lower Quartile"] = V_Q25 + results[%$"Upper Quartile"] = V_Q75 + results[%$"Inter-quartile range"] = V_IQR + results[%$"Median absolute deviation"] = V_MAD + results[%$"Most frequent value"] = V_Mode endif WAVE/T statsLabels = ListToTextWave(PSX_STATS_LABELS, ";") @@ -1710,24 +1862,39 @@ static Function/WAVE PSX_OperationStatsImpl(string graph, string id, WAVE/WAVE s return output End -threadsafe static Function PSX_CalculateRiseTime(WAVE sweepDataOffFilt, WAVE psxEvent, variable kernelAmp, variable lowerThreshold, variable upperThreshold, variable index) +threadsafe static Function PSX_CalculateRiseTimeWrapper(WAVE sweepDataOffFilt, WAVE psxEvent, variable kernelAmp, variable lowerThreshold, variable upperThreshold, variable index) - variable dY, xStart, xEnd, yStart, yEnd, xlt, xupt, lowerLevel, upperLevel, riseTime - variable printDebug - string comboKey + variable error, xStart, xEnd, yStart, yEnd, riseTime + string comboKey // deconvPeak is defined in the deconvoluted wave, // so we can't use %deconvPeak as y-value xStart = psxEvent[index][%$"Onset Time"] + yStart = IsNaN(xStart) ? NaN : sweepDataOffFilt(xStart) - if(IsNaN(xStart)) - return NaN + xEnd = psxEvent[index][%peak_t] + yEnd = psxEvent[index][%peak] + + [riseTime, error] = PSX_CalculateRiseTime(sweepDataOffFilt, xStart, xEnd, yStart, yEnd, lowerThreshold, upperThreshold, kernelAmp) + +#ifdef DEBUGGING_ENABLED + if(error) + comboKey = JWN_GetStringFromWaveNote(psxEvent, PSX_EVENTS_COMBO_KEY_WAVE_NOTE) + + printf "comboKey: %s, x: [%g, %g], y: [%g, %g], index: %d, thresholds: [%g, %g], risetime: %g\r", comboKey, xStart, xEnd, yStart, yEnd, index, lowerThreshold, upperThreshold, risetime endif +#endif // DEBUGGING_ENABLED - yStart = sweepDataOffFilt(xStart) + return riseTime +End - xEnd = psxEvent[index][%peak_t] - yEnd = psxEvent[index][%peak] +threadsafe static Function [variable riseTime, variable error] PSX_CalculateRiseTime(WAVE data, variable xStart, variable xEnd, variable yStart, variable yEnd, variable lowerThreshold, variable upperThreshold, variable kernelAmp) + + variable dY, xlt, xupt, lowerLevel, upperLevel + + if(IsNaN(xStart)) + return [NaN, 1] + endif dY = abs(yStart - yEnd) @@ -1738,34 +1905,26 @@ threadsafe static Function PSX_CalculateRiseTime(WAVE sweepDataOffFilt, WAVE psx xlt = NaN xupt = NaN - FindLevel/R=(xStart, xEnd)/Q sweepDataOffFilt, lowerLevel + FindLevel/R=(xStart, xEnd)/Q data, lowerLevel if(!V_flag) xlt = V_levelX else - printDebug = 1 + error = 1 endif - FindLevel/R=(xStart, xEnd)/Q sweepDataOffFilt, upperLevel + FindLevel/R=(xStart, xEnd)/Q data, upperLevel if(!V_flag) xupt = V_levelX else - printDebug = 1 + error = 1 endif ASSERT_TS(kernelAmp != 0 && IsFinite(kernelAmp), "kernelAmp must be finite and not zero") riseTime = (xlt - xupt) * sign(kernelAmp) * (-1) -#ifdef DEBUGGING_ENABLED - if(printDebug) - comboKey = JWN_GetStringFromWaveNote(psxEvent, PSX_EVENTS_COMBO_KEY_WAVE_NOTE) - - printf "comboKey: %s, x: [%g, %g], y: [%g, %g], index: %d, dY: %g, thresholds: [%g, %g], levels: [%g, %g], risetime: %g, xlt: %g, xupt: %g\r", comboKey, xStart, xEnd, yStart, yEnd, index, dY, lowerThreshold, upperThreshold, lowerLevel, upperLevel, risetime, xlt, xupt - endif -#endif // DEBUGGING_ENABLED - - return riseTime + return [riseTime, error] End threadsafe static Function PSX_CalculateOnsetTime(WAVE sweepDataDiff, WAVE psxEvent, variable kernelAmp, variable diffThreshPerc, variable index) @@ -1874,6 +2033,8 @@ Function/S PSX_FitResultToString(variable fitResult) return UpperCaseFirstChar(GetErrMessage(abs(fitResult))) elseif(fitResult == PSX_DECAY_FIT_ERROR) return "Too large tau" + elseif(fitResult == PSX_DECAY_FIT_INVALID_RANGE_ERROR) + return "Invalid fit range" endif BUG("Unknown fitResult") @@ -1989,7 +2150,7 @@ End static Function PSX_UpdateOffsetInAllEventGraph(string win) string extAllGraph, specialEventPanel - variable i, numEvents, offsetMode, first, last, xOffset, yOffset + variable i, numEvents, offsetMode, first, last, xOffset extAllGraph = PSX_GetAllEventGraph(win) @@ -2026,26 +2187,22 @@ static Function PSX_UpdateOffsetInAllEventGraph(string win) [first, last] = PSX_GetSingleEventRange(psxEvent, sweepDataOffFilt, i) - Duplicate/FREE/R=(first, last) sweepDataOffFilt, singleEventRaw + Duplicate/O/R=(first, last) sweepDataOffFilt, singleEvent switch(offsetMode) case PSX_HORIZ_OFFSET_ONSET: xOffset = IsFinite(psxEvent[i][%$"Onset Time"]) ? (first - psxEvent[i][%$"Onset Time"]) : 0 - yOffset = 0 break case PSX_HORIZ_OFFSET_PEAK: xOffset = first - psxEvent[i][%peak_t] - yOffset = 0 break case PSX_HORIZ_OFFSET_SLEW: xOffset = first - psxEvent[i][%$"Slew Rate Time"] - yOffset = 0 break default: FATAL_ERROR("Invalid offset mode") endswitch - MultiThread singleEvent[] = singleEventRaw[p] - yOffset SetScale/P x, xOffset, DimDelta(singleEvent, ROWS), singleEvent endfor endfor @@ -2210,26 +2367,28 @@ static Function PSX_AdaptColorsInAllEventGraph(string win, [variable forceAverag endif if(averageUpdateRequired || forceAverageUpdate) - PSX_UpdateAverageTraces(win, eventIndexFromTraces, comboIndizesFromTraces, stateNew, indexMapper, comboFolders) + PSX_UpdateAverageTraces(win, comboIndizesFromTraces, stateNew, indexMapper, comboFolders) endif End /// @brief Update the contents of the average waves for the all event graph -static Function PSX_UpdateAverageTraces(string win, WAVE/T eventIndexFromTraces, WAVE/T comboIndizesFromTraces, WAVE stateNew, WAVE indexMapper, WAVE/DF comboFolders) +static Function PSX_UpdateAverageTraces(string win, WAVE/T comboIndizesFromTraces, WAVE stateNew, WAVE indexMapper, WAVE/DF comboFolders) variable i, idx, numEvents, eventState, start, stop variable acceptIndex, rejectIndex, undetIndex, extractStartAbs, extractStopAbs, fitStartAbs - string extAllGraph, name + string extAllGraph, name, path extAllGraph = PSX_GetAllEventGraph(win) - numEvents = DimSize(eventIndexFromTraces, ROWS) + numEvents = DimSize(stateNew, ROWS) Make/WAVE/FREE/N=(numEvents) contAverageAll, contAverageAccept, contAverageReject, contAverageUndet - Make/FREE/D/N=(numEvents) eventStopTime + Make/FREE=1/D/N=(numEvents) eventOnsetTime, eventPeakTime, eventStopTime, eventKernelAmp for(i = 0; i < numEvents; i += 1) - idx = str2num(eventIndexFromTraces[i]) + // i: index of the displayed events across multiple combos + // idx: index in a single combo + idx = indexMapper[i] DFREF comboDFR = comboFolders[str2num(comboIndizesFromTraces[i])] @@ -2239,18 +2398,22 @@ static Function PSX_UpdateAverageTraces(string win, WAVE/T eventIndexFromTraces, WAVE/SDFR=singleEventDFR singleEvent = $name - switch(stateNew[i]) - case PSX_ACCEPT: - contAverageAccept[acceptIndex] = singleEvent + WAVE sweepDataOffFilt = GetPSXSweepDataOffFiltWaveFromDFR(comboDFR) + WAVE psxEvent = GetPSXEventWaveFromDFR(comboDFR) - WAVE sweepDataOffFilt = GetPSXSweepDataOffFiltWaveFromDFR(comboDFR) - WAVE psxEvent = GetPSXEventWaveFromDFR(comboDFR) + // single event waves are zeroed in x-direction to extractStartAbs + [extractStartAbs, extractStopAbs] = PSX_GetSingleEventRange(psxEvent, sweepDataOffFilt, idx) + eventStopTime[i] = extractStopAbs - extractStartAbs + eventOnsetTime[i] = psxEvent[idx][%$"Onset Time"] - extractStartAbs + eventPeakTime[i] = psxEvent[idx][%peak_t] - extractStartAbs - // single event waves are zeroed in x-direction to extractStartAbs - [extractStartAbs, extractStopAbs] = PSX_GetSingleEventRange(psxEvent, sweepDataOffFilt, idx) - eventStopTime[acceptIndex] = extractStopAbs - extractStartAbs + path = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX_KERNEL + eventKernelAmp[i] = JWN_GetNumberFromWaveNote(psxEvent, path + "/amp") - acceptIndex += 1 + switch(stateNew[i]) + case PSX_ACCEPT: + contAverageAccept[acceptIndex] = singleEvent + acceptIndex += 1 break case PSX_REJECT: contAverageReject[rejectIndex] = singleEvent @@ -2270,12 +2433,16 @@ static Function PSX_UpdateAverageTraces(string win, WAVE/T eventIndexFromTraces, DFREF averageDFR = PSX_GetAverageFolder(win) PSX_UpdateAverageWave(contAverageAccept, acceptIndex, averageDFR, PSX_ACCEPT) + PSX_FitAverage(win, averageDFR, eventOnsetTime, eventPeakTime, eventStopTime, eventKernelAmp, stateNew, PSX_ACCEPT) + PSX_UpdateAverageWave(contAverageReject, rejectIndex, averageDFR, PSX_REJECT) + PSX_FitAverage(win, averageDFR, eventOnsetTime, eventPeakTime, eventStopTime, eventKernelAmp, stateNew, PSX_REJECT) + PSX_UpdateAverageWave(contAverageUndet, undetIndex, averageDFR, PSX_UNDET) - PSX_UpdateAverageWave(contAverageAll, numEvents, averageDFR, PSX_ALL) + PSX_FitAverage(win, averageDFR, eventOnsetTime, eventPeakTime, eventStopTime, eventKernelAmp, stateNew, PSX_UNDET) - Redimension/N=(acceptIndex) eventStopTime - PSX_FitAcceptAverage(win, averageDFR, eventStopTime) + PSX_UpdateAverageWave(contAverageAll, numEvents, averageDFR, PSX_ALL) + PSX_FitAverage(win, averageDFR, eventOnsetTime, eventPeakTime, eventStopTime, eventKernelAmp, stateNew, PSX_ALL) End /// @brief Helper function to update the average waves for the all event graph @@ -2299,82 +2466,325 @@ static Function/DF PSX_GetAverageFolder(string win) return PSX_GetWorkingFolder(win) End -static Function PSX_FitAcceptAverage(string win, DFREF averageDFR, WAVE eventStopTime) +Function PSX_UpdateInfoButtonHelp(string win, WAVE/T help, variable state) + + string ctrl + + ctrl = "button_fit_results_" + LowerStr(PSX_StateToString(state)) + + UpdateInfoButtonHelp(win, ctrl, help) +End + +Function PSX_IsFitEnabled(string win, variable state) + + string ctrl + + ctrl = "checkbox_events_fit_" + LowerStr(PSX_StateToString(state)) + + return GetCheckBoxState(win, ctrl) +End + +Function PSX_FilterFitAverageParameters(WAVE input, WAVE newState, variable state, [variable requireFiniteResult]) + + variable numEvents = DimSize(newState, ROWS) + + if(ParamIsDefault(requireFiniteResult)) + requireFiniteResult = 0 + else + requireFiniteResult = !!requireFiniteResult + endif + + Make/FREE/N=(numEvents) result = (newState[p] & state) ? input[p] : NaN + + WAVE/Z resultClean = ZapNaNs(result) + if(WaveExists(resultClean)) + return mean(resultClean) + endif + + ASSERT(!requireFiniteResult, "Could not find any events with finite " + NameOfWave(input)) + + return Inf +End + +/// @brief Compute and store fitted average traces for a given event state +/// +/// Fits average, rise, and decay components of events marked with the specified state. +/// If fitting is disabled via GUI checkbox, outputs are set to NaN. +/// +/// @param win Name of the calling window (used to access GUI settings) +/// @param averageDFR Data folder reference containing the average waves +/// @param eventOnsetTime Wave of event onset times +/// @param eventPeakTime Wave of event peak times +/// @param eventStopTime Wave of event stop times +/// @param eventKernelAmp Wave of event amplitudes (used to infer polarity) +/// @param newState Wave containing per-event classification bit flags +/// @param state Bit flag representing the target event state (e.g., PSX_ACCEPT) +static Function PSX_FitAverage(string win, DFREF averageDFR, WAVE eventOnsetTime, WAVE eventPeakTime, WAVE eventStopTime, WAVE eventKernelAmp, WAVE newState, variable state) string specialEventPanel, str, htmlStr, rawCode, browser, msg, fitFunc - variable err, numAveragePoints, start, stop, meanStopTime + variable err, numAveragePoints, start, stop, meanStopTime, meanOnsetTime, meanPeakTime, meanKernelAmp + variable xStart, xEnd, yStart, yEnd, riselowerThreshold, riseAndDecayUpperThreshold + variable extrema, extrema_t, edge, riseStart, riseStop, decayStart, decayStop, wTau, backwardEdge - WAVE acceptedAverageFit = GetPSXAcceptedAverageFitWaveFromDFR(averageDFR) + WAVE AverageFit = GetPSXAverageFitWaveFromDFR(averageDFR, state) + WAVE RiseAverageFit = GetPSXRiseAverageFitWaveFromDFR(averageDFR, state) + WAVE DecayAverageFit = GetPSXDecayAverageFitWaveFromDFR(averageDFR, state) specialEventPanel = PSX_GetSpecialPanel(win) Make/T/FREE help = {PSX_AVERAGE_FIT_RESULT_DEFAULT_HELP} - UpdateInfoButtonHelp(specialEventPanel, "button_fit_results", help) + PSX_UpdateInfoButtonHelp(specialEventPanel, help, state) - if(!GetCheckBoxState(specialEventPanel, "checkbox_average_events_fit")) - FastOp acceptedAverageFit = (NaN) + if(!PSX_IsFitEnabled(specialEventPanel, state)) + FastOp AverageFit = (NaN) + FastOp RiseAverageFit = (NaN) + FastOp DecayAverageFit = (NaN) return NaN endif - WAVE average = GetPSXAverageWave(averageDFR, PSX_ACCEPT) + WAVE average = GetPSXAverageWave(averageDFR, state) numAveragePoints = DimSize(average, ROWS) if(numAveragePoints == 0 || !HasOneValidEntry(average)) - FastOp acceptedAverageFit = (NaN) + FastOp AverageFit = (NaN) + FastOp RiseAverageFit = (NaN) + FastOp DecayAverageFit = (NaN) return NaN endif - Redimension/N=(numAveragePoints) acceptedAverageFit - FastOp acceptedAverageFit = (NaN) - CopyScales average, acceptedAverageFit - - WAVE/Z eventStopTimeClean = ZapNaNs(eventStopTime) - if(WaveExists(eventStopTimeClean)) - meanStopTime = mean(eventStopTime) + Redimension/N=(numAveragePoints) AverageFit, RiseAverageFit, DecayAverageFit + FastOp AverageFit = (NaN) + FastOp RiseAverageFit = (NaN) + FastOp DecayAverageFit = (NaN) + CopyScales average, AverageFit, RiseAverageFit, DecayAverageFit + + meanOnsetTime = PSX_FilterFitAverageParameters(eventOnsetTime, newState, state) + // meanStopTime = PSX_FindDominantStopTimePeak(eventStopTime, newState, state, eventKernelAmp, binwidth = 0.1, showplot = 1) + meanStopTime = PSX_FindDominantStopTimePeak(eventStopTime, newState, state, eventKernelAmp) + meanKernelAmp = PSX_FilterFitAverageParameters(eventKernelAmp, newState, state, requireFiniteResult = 1) + meanPeakTime = PSX_FilterFitAverageParameters(eventPeakTime, newState, state, requireFiniteResult = 1) + + WaveStats/M=1/Q average + + if(meanKernelAmp > 0) + extrema = V_max + extrema_t = V_maxLoc + edge = FINDLEVEL_EDGE_INCREASING + backwardEdge = FINDLEVEL_EDGE_DECREASING + elseif(meanKernelAmp < 0) + extrema = V_min + extrema_t = V_minLoc + edge = FINDLEVEL_EDGE_DECREASING + backwardEdge = FINDLEVEL_EDGE_INCREASING else - meanStopTime = Inf + FATAL_ERROR("Invalid kernel amp") endif - start = 0 - stop = min(IndexToScale(average, DimSize(average, ROWS) - 1, ROWS), meanStopTime) + riselowerThreshold = GetSetVariable(specialEventPanel, "setvar_fit_start_amplitude") * PERCENT_TO_ONE + riseAndDecayUpperThreshold = 0.9 + // FindLevel/EDGE=(edge)/Q average, (riselowerThreshold * extrema) + variable onsetX = PSX_CalculateOnsetTimeFromAvg(average, meanKernelAmp, meanOnsetTime, meanPeakTime) + assert(!isNaN(onsetX), "onset time calculation must cannot be Nan") + variable level = (riseLowerThreshold * (extrema - average(onsetX))) + average(onsetX) + FindLevel/EDGE=(edge)/R=(extrema_t, onsetX)/Q average, level + riseStart = V_LevelX + FindLevel/EDGE=(edge)/Q average, (riseAndDecayUpperThreshold * extrema) + riseStop = V_levelX + + FindLevel/EDGE=(backwardEdge)/R=(Inf, 0)/Q average, (riseAndDecayUpperThreshold * extrema) + decayStart = V_levelX + + decayStop = min(IndexToScale(average, DimSize(average, ROWS) - 1, ROWS), meanStopTime) AssertOnAndClearRTError() - fitFunc = GetPopupMenuString(specialEventPanel, "popupmenu_accept_fit_function") - strswitch(fitFunc) - case "dblexp_peak": - Make/FREE/D/N=5 coefWave - Make/FREE/T coeffNames = {"y0", "A", "tau1", "tau2", "X0"} - CurveFit/M=0/Q/N=2 dblexp_peak, kwCWave=coefWave, average(start, stop)/D=acceptedAverageFit; err = GetRTError(1) - break - case "dblexp_XOffset": - Make/FREE/D/N=5 coefWave - Make/FREE/T coeffNames = {"y0", "A1", "tau1", "A2", "tau2"} - CurveFit/M=0/Q/N=2 dblexp_XOffset, kwCWave=coefWave, average(start, stop)/D=acceptedAverageFit; err = GetRTError(1) - break - default: - FATAL_ERROR("Unknown fit function") - endswitch - sprintf msg, "Fit in the range [%g, %g] finished with %d (%s)\r", start, stop, err, GetErrMessage(err) + Make/FREE/T/N=(11, 2) InputAvg, InputRise, InputDecay + + Make/FREE/D/N=5 coefWave + Make/FREE/T coeffNames = {"y0", "A", "tau1", "tau2", "X0"} + CurveFit/M=0/Q/N=2 dblexp_peak, kwCWave=coefWave, average(riseStart, decayStop)/D=AverageFit; err = GetRTError(1) + + InputAvg[0][0, 1] = {{"Function"}, {"dblexp_peak"}} + InputAvg[1][0, 1] = {{"ChiSq"}, {num2strHighPrec(V_chiSq)}} + InputAvg[2, 6][0] = coeffNames[p - 2] + InputAvg[2, 6][1] = num2strHighPrec(coefWave[p - 2]) + InputAvg[7][0, 1] = {{"Maximum"}, {num2strHighPrec(WaveMax(AverageFit))}} + InputAvg[8][0, 1] = {{"State source"}, {PSX_GetStateTypeFromSpecialPanel(win)}} + InputAvg[9][0, 1] = {{"Current combo"}, {ToTrueFalse(PSX_GetrestrictEventsToCurrentCombo(win))}} + + Make/FREE/D/N=5 coefWave + Make/FREE/T coeffNames = {"y0", "A1", "tau1", "A2", "tau2"} + CurveFit/M=0/Q/N=1 dblexp_XOffset, kwCWave=coefWave, average(riseStart, riseStop)/D=RiseAverageFit; err = GetRTError(1) + + wTau = ((coefWave[1] * coefWave[2] + coefWave[3] * coefWave[4]) / (coefWave[1] + coefWave[3])) + + InputRise[0][0, 1] = {{"Function"}, {"dblexp_XOffset rise"}} + InputRise[1][0, 1] = {{"ChiSq"}, {num2strHighPrec(V_chiSq)}} + InputRise[2, 6][0] = coeffNames[p - 2] + InputRise[2, 6][1] = num2strHighPrec(coefWave[p - 2]) + InputRise[7][0, 1] = {{"Maximum"}, {num2strHighPrec(WaveMax(RiseAverageFit))}} + InputRise[8][0, 1] = {{"weighted tau"}, {num2strHighPrec(wTau)}} + InputRise[9][0, 1] = {{"State source"}, {PSX_GetStateTypeFromSpecialPanel(win)}} + InputRise[10][0, 1] = {{"Current combo"}, {ToTrueFalse(PSX_GetrestrictEventsToCurrentCombo(win))}} + + Make/FREE/D/N=5 coefWave + Make/FREE/T coeffNames = {"y0", "A1", "tau1", "A2", "tau2"} + CurveFit/M=0/Q/N=1 dblexp_XOffset, kwCWave=coefWave, average(decayStart, decayStop)/D=DecayAverageFit; err = GetRTError(1) + + wTau = ((coefWave[1] * coefWave[2] + coefWave[3] * coefWave[4]) / (coefWave[1] + coefWave[3])) + + InputDecay[0][0, 1] = {{"Function"}, {"dblexp_XOffset decay"}} + InputDecay[1][0, 1] = {{"ChiSq"}, {num2strHighPrec(V_chiSq)}} + InputDecay[2, 6][0] = coeffNames[p - 2] + InputDecay[2, 6][1] = num2strHighPrec(coefWave[p - 2]) + InputDecay[7][0, 1] = {{"Maximum"}, {num2strHighPrec(WaveMax(DecayAverageFit))}} + InputDecay[8][0, 1] = {{"weighted tau"}, {num2strHighPrec(wTau)}} + InputDecay[9][0, 1] = {{"State source"}, {PSX_GetStateTypeFromSpecialPanel(win)}} + InputDecay[10][0, 1] = {{"Current combo"}, {ToTrueFalse(PSX_GetrestrictEventsToCurrentCombo(win))}} + + sprintf msg, "Fit in the range [%g, %g] finished with %d (%s)\r", riseStart, decayStop, err, GetErrMessage(err) DEBUGPRINT(msg) if(err) return NaN endif - Make/FREE/T/N=(9, 2) input + Concatenate/NP/FREE {InputAvg, InputRise, InputDecay}, input - input[0][0, 1] = {{"Function"}, {fitFunc}} - input[1][0, 1] = {{"ChiSq"}, {num2strHighPrec(V_chiSq)}} - input[2, 6][0] = coeffNames[p - 2] - input[2, 6][1] = num2strHighPrec(coefWave[p - 2]) - input[7][0, 1] = {{"State source"}, {PSX_GetStateTypeFromSpecialPanel(win)}} - input[8][0, 1] = {{"Current combo"}, {ToTrueFalse(PSX_GetrestrictEventsToCurrentCombo(win))}} - - UpdateInfoButtonHelp(specialEventPanel, "button_fit_results", input) + PSX_UpdateInfoButtonHelp(specialEventPanel, input, state) browser = SFH_GetBrowserForFormulaGraph(win) PSX_StoreIntoResultsWave(browser, SFH_RESULT_TYPE_PSX_MISC, input, "accepted average fit results") + + WAVE AverageFitResults = GetPSXAverageFitResultsWaveFromDFR(averageDFR, state) + Redimension/N=(0, -1) AverageFitResults + + Concatenate/NP=(ROWS) {InputAvg, InputRise, InputDecay}, AverageFitResults +End + +//// @brief Return most frequent stop time peak, ignoring edge bins and bins < decayTau. Falls back to 2×decayTau. +/// @param stopTimes Wave of event stop times (with kernel decayTau in wave note) +/// @param newState Bitmask wave of per-event state +/// @param state Target state value (e.g., PSX_ACCEPT) +/// @param psxKernel kernel +/// @param [binWidth] Histogram bin width (default = 0.1) +/// @param [showPlot] Whether to display the histogram (default = 0) +Function PSX_FindDominantStopTimePeak(WAVE stopTimes, WAVE newState, variable state, WAVE psxKernel, [variable binWidth, variable showPlot]) + + variable i + variable n = DimSize(stopTimes, ROWS) + variable count = 0 + + if(ParamIsDefault(binWidth)) + binWidth = 0.1 + endif + if(ParamIsDefault(showPlot)) + showPlot = 0 + endif + + // Step 1: Filter stopTimes by state + Make/FREE/D/N=(n) tmp = NaN + for(i = 0; i < n; i += 1) + if((newState[i] & state) && IsFinite(stopTimes[i])) + tmp[count] = stopTimes[i] + count += 1 + endif + endfor + + if(count == 0) + if(showPlot) + Print "No valid stop times for selected state." + endif + return NaN + endif + + Redimension/N=(count) tmp + + // Step 2: Bin setup + variable minVal = WaveMin(tmp) + variable maxVal = WaveMax(tmp) + + if(minVal == maxVal) + if(showPlot) + Print "All stop times identical. Using ", minVal + endif + return minVal + endif + + variable nBins = ceil((maxVal - minVal) / binWidth) + if(nBins < 5) + if(showPlot) + Print "Too few bins. Returning NaN." + endif + return NaN + endif + + // Step 3: Create histogram + KillWaves/Z root:PSX_hist, root:PSX_histX + Make/O/D/N=(nBins) root:PSX_hist, root:PSX_histX + WAVE hist = root:PSX_hist + WAVE histX = root:PSX_histX + + Histogram/B={minVal, binWidth, nBins}/DEST=hist tmp + histX[] = minVal + binWidth * (p + 0.5) + Smooth 3, hist + + // Step 4: Optional plot + if(showPlot) + // Display/K=1/N=PSX_StopTimeHistogram + + if(!WindowExists("PSX_StopTimeHistogram")) // 1 = graph window + Display/N=PSX_StopTimeHistogram + else + DoWindow/F PSX_StopTimeHistogram // Bring to front if already exists + // AppendToGraph/W=PSX_StopTimeHistogram stopTimes vs stopTimeHist + endif + + AppendToGraph/W=PSX_StopTimeHistogram hist vs histX + ModifyGraph mode=4, lsize=2, rgb=(0, 0, 65535) + ModifyGraph grid=1, mirror=1, axThick=1.5 + Label bottom, "Stop Time" + Label left, "Event Count" + SetAxis/A + endif + + // Step 5: Try to find decayTau + variable decayTau = GetKernelDecayTau() + variable validTau = IsFinite(decayTau) + + // Step 6: Find max internal peak above decayTau + variable peakMaxVal = -Inf + variable maxIdx = -1 + for(i = 1; i < (nBins - 1); i += 1) + if(validTau && histX[i] < (2 * decayTau)) + continue + endif + if(hist[i] > peakMaxVal) + peakMaxVal = hist[i] + maxIdx = i + endif + endfor + + // Step 7: Return result + if(maxIdx >= 0) + if(showPlot) + Print "Dominant peak at ", histX[maxIdx], " (passes decayTau filter)" + endif + return histX[maxIdx] + endif + + // Step 8: Fallback + if(validTau) + if(showPlot) + Print "No peak ≥ decayTau found. Fallback to 2 × decayTau = ", 2 * decayTau + endif + return 2 * decayTau + endif + + if(showPlot) + Print "No peak found. DecayTau not available. Returning NaN." + endif + return NaN End static Function PSX_StoreIntoResultsWave(string browser, variable resultType, WAVE data, string name) @@ -2747,7 +3157,7 @@ End static Function PSX_AppendAverageTraces(string extAllGraph, DFREF averageDFR, string traceSuffix, variable idx, string comboKey, variable comboIndex, WAVE traceUserDataKeys, WAVE states, WAVE acceptColors, WAVE rejectColors, WAVE undetColors) variable state - string trace + string trace, stateAsString, nameTemplate for(state : states) @@ -2766,18 +3176,53 @@ static Function PSX_AppendAverageTraces(string extAllGraph, DFREF averageDFR, st TUD_SetUserDataFromWaves(extAllGraph, trace, \ traceUserDataKeys, \ {"NaN", num2str(state), num2str(state), "0", PSX_TUD_TYPE_AVERAGE, comboKey, "NaN", num2str(comboIndex)}) + + idx += 1 + + stateAsString = PSX_StateToString(state) + + // averageFit + WAVE averageFit = GetPSXAverageFitWaveFromDFR(averageDFR, state) + + nameTemplate = stateAsString + "_AverageFit" + trace = PSX_GetAverageTraceName(idx, nameTemplate, comboIndex, traceSuffix) + + AppendToGraph/W=$extAllGraph averageFit/TN=$trace + + TUD_SetUserDataFromWaves(extAllGraph, trace, \ + traceUserDataKeys, \ + {"NaN", num2str(state), num2str(state), "0", PSX_TUD_TYPE_AVERAGE, comboKey, "NaN", num2str(comboIndex)}) + idx += 1 - endfor - WAVE acceptedAverageFit = GetPSXAcceptedAverageFitWaveFromDFR(averageDFR) + // riseAverageFit + WAVE riseAverageFit = GetPSXRiseAverageFitWaveFromDFR(averageDFR, state) - trace = PSX_GetAverageTraceName(idx, "acceptAverageFit", comboIndex, traceSuffix) - idx += 1 + nameTemplate = stateAsString + "_RiseAverageFit" + trace = PSX_GetAverageTraceName(idx, nameTemplate, comboIndex, traceSuffix) - AppendToGraph/W=$extAllGraph acceptedAverageFit/TN=$trace - TUD_SetUserDataFromWaves(extAllGraph, trace, \ - traceUserDataKeys, \ - {"NaN", num2str(PSX_ACCEPT), num2str(PSX_ACCEPT), "0", PSX_TUD_TYPE_AVERAGE, comboKey, "NaN", num2str(comboIndex)}) + AppendToGraph/W=$extAllGraph riseAverageFit/TN=$trace + + TUD_SetUserDataFromWaves(extAllGraph, trace, \ + traceUserDataKeys, \ + {"NaN", num2str(state), num2str(state), "0", PSX_TUD_TYPE_AVERAGE, comboKey, "NaN", num2str(comboIndex)}) + + idx += 1 + + // decayAverageFit + WAVE decayAverageFit = GetPSXDecayAverageFitWaveFromDFR(averageDFR, state) + + nameTemplate = stateAsString + "_DecayAverageFit" + trace = PSX_GetAverageTraceName(idx, nameTemplate, comboIndex, traceSuffix) + + AppendToGraph/W=$extAllGraph decayAverageFit/TN=$trace + + TUD_SetUserDataFromWaves(extAllGraph, trace, \ + traceUserDataKeys, \ + {"NaN", num2str(state), num2str(state), "0", PSX_TUD_TYPE_AVERAGE, comboKey, "NaN", num2str(comboIndex)}) + + idx += 1 + endfor return idx End @@ -2833,9 +3278,9 @@ static Function PSX_UpdateSingleEventTextbox(string win, [variable eventIndex]) DFREF comboDFR = PSX_GetCurrentComboFolder(win) WAVE psxEvent = GetPSXEventWaveFromDFR(comboDFR) - Make/FREE/T/N=(8, 2) input + Make/FREE/T/N=(10, 2) input - input[0][0] = {"Event State:", "Fit State:", "Fit Result:", "Event:", "Deconv Peak:", "Peak:", "Baseline:", "IeI:", "Amp (rel.):", "Tau:", "Onset time:", "Rise time:"} + input[0][0] = {"Event State:", "Fit State:", "Fit Result:", "Event:", "Deconv Peak:", "Peak:", "Baseline:", "IeI:", "Amp (rel.):", "Slow tau:", "Fast tau:", "Weighted tau:", "Onset time:", "Rise time:"} input[0][1] = {PSX_StateToString(psxEvent[eventIndex][%$"Event manual QC call"]), \ PSX_StateToString(psxEvent[eventIndex][%$"Fit manual QC call"]), \ PSX_FitResultToString(psxEvent[eventIndex][%$"Fit Result"]), \ @@ -2845,7 +3290,9 @@ static Function PSX_UpdateSingleEventTextbox(string win, [variable eventIndex]) num2str(psxEvent[eventIndex][%baseline_t], "%8.02f") + " [ms]", \ num2str(psxEvent[eventIndex][%iei], "%8.02f") + " [ms]", \ num2str(psxEvent[eventIndex][%amplitude], "%8.02f") + " [" + yUnit + "]", \ - num2str(psxEvent[eventIndex][%tau], "%8.02f") + " [ms]", \ + num2str(psxEvent[eventIndex][%slowTau], "%8.02f") + " [ms]", \ + num2str(psxEvent[eventIndex][%fastTau], "%8.02f") + " [ms]", \ + num2str(psxEvent[eventIndex][%weightedTau], "%8.02f") + " [ms]", \ num2str(psxEvent[eventIndex][%$"Onset Time"], "%8.02f") + " [ms]", \ num2str(psxEvent[eventIndex][%$"Rise Time"], "%8.02f") + " [ms]"} @@ -3546,9 +3993,16 @@ static Function PSX_StoreGuiState(string win, string browser) JSON_SetVariable(jsonID, "/specialEventPanel/popup_block", GetPopupMenuIndex(specialEventPanel, "popup_block")) JSON_SetVariable(jsonID, "/specialEventPanel/setvar_event_block_size", GetSetVariable(specialEventPanel, "setvar_event_block_size")) + JSON_SetVariable(jsonID, "/specialEventPanel/setvar_fit_start_amplitude", GetSetVariable(specialEventPanel, "setvar_fit_start_amplitude")) mainWindow = GetMainWindow(win) - JSON_SetVariable(jsonID, "/mainPanel/checkbox_suppress_update", GetCheckBoxState(mainWindow, "checkbox_suppress_update")) + + WAVE/T checkboxes = PSX_GetSpecialEventPanelCheckboxes(mainWindow) + + for(ctrl : checkboxes) + JSON_SetVariable(jsonID, "/mainPanel/" + ctrl, GetCheckBoxState(mainWindow, ctrl)) + endfor + JSON_SetVariable(jsonID, "/mainPanel/listbox_select_combo", GetListBoxSelRow(mainWindow, "listbox_select_combo")) WAVE axesProps = GetAxesProperties(extAllGraph) @@ -3590,22 +4044,23 @@ static Function PSX_RestoreGuiState(string win) SetCheckBoxState(specialEventPanel, ctrl, JSON_GetVariable(jsonID, "/specialEventPanel/" + ctrl)) endfor - // first block size, as that recalculates the number of blocks - PGC_SetAndActivateControl(specialEventPanel, "setvar_event_block_size", val = JSON_GetVariable(jsonID, "/specialEventPanel/setvar_event_block_size")) + SetSetVariable(specialEventPanel, "setvar_event_block_size", JSON_GetVariable(jsonID, "/specialEventPanel/setvar_event_block_size")) + SetSetVariable(specialEventPanel, "setvar_fit_start_amplitude", JSON_GetVariable(jsonID, "/specialEventPanel/setvar_fit_start_amplitude")) - selectedBlock = JSON_GetVariable(jsonID, "/specialEventPanel/popup_block") - allBlocks = PSX_GetAllEventBlockNumbers(specialEventPanel) - lastBlock = NumberFromList(ItemsInList(allBlocks) - 1, allBlocks) - PGC_SetAndActivateControl(specialEventPanel, "popup_block", val = limit(selectedBlock, 0, lastBlock)) + SetPopupMenuIndex(specialEventPanel, "popup_block", JSON_GetVariable(jsonID, "/specialEventPanel/popup_block")) WAVE/T popups = PSX_GetSpecialEventPanelPopups(specialEventPanel) for(popMenu : popups) - PGC_SetAndActivateControl(specialEventPanel, popMenu, val = JSON_GetVariable(jsonID, "/specialEventPanel/" + popMenu)) + SetPopupMenuIndex(specialEventPanel, popMenu, JSON_GetVariable(jsonID, "/specialEventPanel/" + popMenu)) endfor mainWindow = GetMainWindow(win) - SetCheckBoxState(mainWindow, "checkbox_suppress_update", JSON_GetVariable(jsonID, "/mainPanel/checkbox_suppress_update")) + WAVE/T controls = PSX_GetSpecialEventPanelCheckboxes(mainWindow) + + for(ctrl : controls) + SetCheckBoxState(mainWindow, ctrl, JSON_GetVariable(jsonID, "/mainPanel/" + ctrl)) + endfor lastActiveCombo = JSON_GetVariable(jsonID, "/mainPanel/listbox_select_combo") @@ -3757,12 +4212,26 @@ static Function PSX_MoveWavesToDataFolders(DFREF workDFR, WAVE/Z/WAVE results, v Duplicate peakY, dfr:peakYAtFilt/WAVE=peakYAtFilt peakYAtFilt[] = sweepDataOffFilt(peakX[p]) - Make/T/N=(numEvents, 2) dfr:eventLocationLabels/WAVE=eventLocationLabels - SetDimLabel COLS, 1, $"Tick Type", eventLocationLabels - eventLocationLabels[][1] = "Major" + WAVE eventLabelsDeconv = GetPSXEventLabelsAsFree(numEvents, PSX_FREE_AXIS_DECONV) + MoveWave eventLabelsDeconv, dfr + + WAVE eventTicksDeconv = GetPSXEventTicksAsFree(numEvents, PSX_FREE_AXIS_DECONV) + MoveWave eventTicksDeconv, dfr + eventTicksDeconv[] = psxEvent[p][%deconvPeak_t] + + WAVE eventLabelsPeak = GetPSXEventLabelsAsFree(numEvents, PSX_FREE_AXIS_PEAK) + MoveWave eventLabelsPeak, dfr - Make/D/N=(numEvents) dfr:eventLocationTicks/WAVE=eventLocationTicks - eventLocationTicks[] = peakX[p] + WAVE eventTicksPeak = GetPSXEventTicksAsFree(numEvents, PSX_FREE_AXIS_PEAK) + MoveWave eventTicksPeak, dfr + eventTicksPeak[] = psxEvent[p][%peak_t] + + WAVE eventLabelsBaseline = GetPSXEventLabelsAsFree(numEvents, PSX_FREE_AXIS_BASELINE) + MoveWave eventLabelsBaseline, dfr + + WAVE eventTicksBaseline = GetPSXEventTicksAsFree(numEvents, PSX_FREE_AXIS_BASELINE) + MoveWave eventTicksBaseline, dfr + eventTicksBaseline[] = psxEvent[p][%baseline_t] PSX_CreateSingleEventWaves(dfr, psxEvent, sweepDataOffFilt) @@ -3909,19 +4378,12 @@ static Function PSX_CreatePSXGraphAndSubwindows(string win, string graph, STRUCT PSX_MarkGraphForPSX(win) - WAVE eventLocationLabels = GetPSXEventLocationLabels(comboDFR) - WAVE eventLocationTicks = GetPSXEventLocationTicks(comboDFR) - WAVE eventColors = GetPSXEventColorsWaveFromDFR(comboDFR) - WAVE eventMarker = GetPSXEventMarkerWaveFromDFR(comboDFR) + WAVE eventColors = GetPSXEventColorsWaveFromDFR(comboDFR) + WAVE eventMarker = GetPSXEventMarkerWaveFromDFR(comboDFR) - NewFreeAxis/W=$win/O/T eventLocAxis - ModifyFreeAxis/W=$win/Z eventLocAxis, master=bottom - ModifyGraph/W=$win grid(eventLocAxis)=1 - ModifyGraph/W=$win tick(eventLocAxis)=3 - ModifyGraph/W=$win lblPos(eventLocAxis)=43 - ModifyGraph/W=$win noLabel(eventLocAxis)=2 - ModifyGraph/W=$win freePos(eventLocAxis)={0, kwFraction} - Label/W=$win eventLocAxis, "\\u#2" + PSX_AppendVisualizationAxis(win, PSX_FREE_AXIS_DECONV) + PSX_AppendVisualizationAxis(win, PSX_FREE_AXIS_PEAK) + PSX_AppendVisualizationAxis(win, PSX_FREE_AXIS_BASELINE) ModifyGraph/W=$win zColor(peakYAtFilt)={eventColors, *, *, directRGB, 0} ModifyGraph/W=$win mode(peakYAtFilt)=3 @@ -3952,6 +4414,47 @@ static Function PSX_CreatePSXGraphAndSubwindows(string win, string graph, STRUCT PS_InitCoordinates(JSONid, mainWin, recursive = 1) End +static Function/WAVE PSX_GetColorsForFreeAxis(string axis) + + string list + + strswitch(axis) + case PSX_FREE_AXIS_BASELINE: + list = PSX_COLORS_BASELINE + break + case PSX_FREE_AXIS_DECONV: + list = PSX_COLORS_DECONV + break + case PSX_FREE_AXIS_PEAK: + list = PSX_COLORS_PEAK + break + default: + FATAL_ERROR("Invalid axis type") + endswitch + + WAVE wv = ListToNumericWave(list, ",") + + return wv +End + +/// @brief Appends a free axis overlayed with the bottom axis name +static Function PSX_AppendVisualizationAxis(string win, string axis) + + NewFreeAxis/W=$win/O/T $axis + ModifyFreeAxis/W=$win/Z $axis, master=bottom + ModifyGraph/W=$win grid($axis)=1 + ModifyGraph/W=$win tick($axis)=3 + ModifyGraph/W=$win lblPos($axis)=43 + ModifyGraph/W=$win noLabel($axis)=2 + ModifyGraph/W=$win freePos($axis)={0, kwFraction} + ModifyGraph/W=$win nticks($axis)=0 + + WAVE RGBcolors = PSX_GetColorsForFreeAxis(axis) + ModifyGraph/W=$win gridRGB($axis)=(RGBcolors[0], RGBcolors[1], RGBcolors[2]) + + Label/W=$win $axis, "\\u#2" +End + /// @brief Mark `win` as being an psx graph static Function PSX_MarkGraphForPSX(string win) @@ -3959,9 +4462,33 @@ static Function PSX_MarkGraphForPSX(string win) End /// @brief Apply plot properties which have to be reapplied on every combo index change -static Function PSX_ApplySpecialPlotProperties(string win, WAVE eventLocationTicks, WAVE eventLocationLabels) +static Function PSX_ApplySpecialPlotProperties(string win) - ModifyGraph/W=$win userticks(eventLocAxis)={eventLocationTicks, eventLocationLabels} + DFREF comboDFR = PSX_GetCurrentComboFolder(win) + + PSX_ApplySpecialPlotPropertiesImpl(win, comboDFR, PSX_FREE_AXIS_DECONV, "checkbox_show_deconv_lines") + PSX_ApplySpecialPlotPropertiesImpl(win, comboDFR, PSX_FREE_AXIS_BASELINE, "checkbox_show_baseline_lines") + PSX_ApplySpecialPlotPropertiesImpl(win, comboDFR, PSX_FREE_AXIS_PEAK, "checkbox_show_peak_lines") +End + +static Function PSX_ApplySpecialPlotPropertiesImpl(string win, DFREF comboDFR, string name, string checkBoxControl) + + variable enabled + string panel, psxGraph + + panel = GetMainWindow(win) + psxGraph = PSX_GetPSXGraph(win) + + enabled = GetCheckboxState(panel, checkBoxControl) + + WAVE eventLabels = GetPSXEventLabelsFromDFR(comboDFR, name) + WAVE eventTicks = GetPSXEventTicksFromDFR(comboDFR, name) + + if(enabled) + ModifyGraph/W=$psxGraph userticks($name)={eventTicks, eventLabels} + else + ModifyGraph/W=$psxGraph userticks($name)=0, nticks($name)=0 + endif End /// @brief Read the user JWN from results and create a legend from all operation parameters @@ -4067,7 +4594,7 @@ End static Function [variable eventIndex, variable waveIndex, variable comboIndex] PSX_GetEventIndexAndComboIndex(string win, [variable direction]) string psxGraph, info, trace, postProc - variable idx, yPointNumber, numEntries + variable idx, yPointNumber, yNumRows, xNumRows psxGraph = PSX_GetPSXGraph(win) @@ -4122,16 +4649,28 @@ static Function [variable eventIndex, variable waveIndex, variable comboIndex] P WAVE/T comboKeys = JWN_GetTextWaveFromWaveNote(yWave, SF_META_USER_GROUP + PSX_JWN_COMBO_KEYS_NAME) - numEntries = DimSize(yWave, ROWS) - ASSERT(numEntries == DimSize(xWave, ROWS), "Unmatching wave sizes") + yNumRows = DimSize(yWave, ROWS) + xNumRows = DimSize(xWave, ROWS) + + if(xNumRows == yNumRows) + // x wave is from the psx yWave + WAVE refxWave = xWave + else + // yWave vs something else + WAVE refxWave = JWN_GetNumericWaveFromWaveNote(yWave, SF_META_XVALUES) + xNumRows = DimSize(refxWave, ROWS) + endif + WaveClear xWave + + ASSERT(xNumRows == yNumRows, "Unmatching wave sizes") - yPointNumber = limit(yPointNumber + direction, 0, numEntries - 1) + yPointNumber = limit(yPointNumber + direction, 0, yNumRows - 1) comboIndex = PSX_GetComboIndexForComboKey(win, comboKeys[yPointNumber]) strswitch(postProc) case "nothing": // fallthrough case "log10": - eventIndex = xWave[yPointNumber] + eventIndex = refxWave[yPointNumber] break case "nonfinite": eventIndex = yWave[yPointNumber] @@ -4231,7 +4770,10 @@ Function PSX_PlotInteractionHook(STRUCT WMWinHookStruct &s) direction = PSX_GetDirectionFromKeyCode(psxGraph, s.keyCode) [eventIndex, waveIndex, comboIndex] = PSX_GetEventIndexAndComboIndex(win, direction = direction) - ASSERT(IsFinite(eventIndex) && IsFinite(waveIndex) && IsFinite(comboIndex), "Invalid event index") + if(IsNaN(eventIndex) || IsNaN(waveIndex) || IsNaN(comboIndex)) + BUG("Could not get new eventIndex after direction application") + break + endif if(PSX_GetCurrentComboIndex(win) != comboIndex) mainWindow = GetMainWindow(win) @@ -4557,11 +5099,25 @@ Function PSX_PostPlot(string win) PSX_UpdateAllEventGraph(win, forceAverageUpdate = 1, forceSingleEventUpdate = 1, forceBlockIndexUpdate = 1, forceSingleEventOffsetUpdate = 1) - DFREF comboDFR = PSX_GetCurrentComboFolder(win) - WAVE eventLocationLabels = GetPSXEventLocationLabels(comboDFR) - WAVE eventLocationTicks = GetPSXEventLocationTicks(comboDFR) + PSX_SetCheckboxTitleColors(win) - PSX_ApplySpecialPlotProperties(win, eventLocationTicks, eventLocationLabels) + PSX_ApplySpecialPlotProperties(win) +End + +static Function PSX_SetCheckboxTitleColors(string win) + + string str, mainWindow + + mainWindow = GetMainWindow(win) + + sprintf str, "Show Deconv \\[0\\K(%s)lines\\]0", PSX_COLORS_DECONV + SetControlTitle(mainWindow, "checkbox_show_deconv_lines", str) + + sprintf str, "Show Peak \\[0\\K(%s)lines\\]0", PSX_COLORS_PEAK + SetControlTitle(mainWindow, "checkbox_show_peak_lines", str) + + sprintf str, "Show Baseline \\[0\\K(%s)lines\\]0", PSX_COLORS_BASELINE + SetControlTitle(mainWindow, "checkbox_show_baseline_lines", str) End static Function PSX_OperationSetDimensionLabels(WAVE/WAVE output, variable numCombos, WAVE/T labels, WAVE/T labelsTemplate) @@ -4641,12 +5197,13 @@ End // Output[1] = sweepDataOffFilt(1) // ... // -// psx(id, [psxKernel(...), numSDs, sweepFilterLow, sweepFilterHigh, maxTauFactor, psxRiseTime(...), psxDeconvFilter(...)]) +// psx(id, [psxKernel(...), numSDs, psxSweepBPFilter(...), maxTauFactor, psxRiseTime(...), psxDeconvBPFilter(...)]) Function/WAVE PSX_Operation(STRUCT SF_ExecutionData &exd) - variable numberOfSDs, sweepFilterLow, sweepFilterHigh, parameterJsonID, numCombos, i, addedData, kernelAmp - variable maxTauFactor, peakThresh, idx, success, kernelRiseTau, kernelDecayTau - string parameterPath, id, psxParameters, dataUnit, path + variable numberOfSDs, parameterJsonID, numCombos, i, readIndex, writeIndex, addedData, kernelAmp + variable maxTauFactor, peakThresh, success, kernelRiseTau, kernelDecayTau, constFilterOrderSweep, constFilterOrderDeconv + variable minFilterOrderSweep, minFilterOrderDeconv + string parameterPath, id, psxParameters, dataUnit, path, msg SFH_CheckArgumentCount(exd, SF_OP_PSX, 1, maxArgs = 8) @@ -4654,69 +5211,116 @@ Function/WAVE PSX_Operation(STRUCT SF_ExecutionData &exd) WAVE/WAVE psxKernelDataset = SFH_GetArgumentAsWave(exd, SF_OP_PSX, 1, defOp = "psxKernel()") + numberOfSDs = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX, 2, defValue = PSX_NUMBER_OF_SDS_DEFAULT, checkFunc = IsStrictlyPositiveAndFinite) + WAVE sweepFilter = SFH_GetArgumentAsWave(exd, SF_OP_PSX, 3, defOp = "psxSweepBPFilter()", singleResult = 1) + maxTauFactor = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX, 5, defValue = PSX_DEFAULT_MAX_TAU_FACTOR, checkFunc = IsStrictlyPositiveAndFinite) + WAVE riseTime = SFH_GetArgumentAsWave(exd, SF_OP_PSX, 6, defOp = "psxRiseTime()", expectedMinorType = IGOR_TYPE_64BIT_FLOAT, singleResult = 1) + ASSERT(IsNumericWave(riseTime), "Invalid return from psxRiseTime") + WAVE deconvFilter = SFH_GetArgumentAsWave(exd, SF_OP_PSX, 7, defOp = "psxDeconvBPFilter()", singleResult = 1) + + path = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX_KERNEL + kernelRiseTau = JWN_GetNumberFromWaveNote(psxKernelDataset, path + "/riseTau") + ASSERT(IsFinite(kernelRiseTau), "riseTau must be finite") + kernelDecayTau = JWN_GetNumberFromWaveNote(psxKernelDataset, path + "/decayTau") + ASSERT(IsFinite(kernelDecayTau), "decayTau must be finite") + + PSX_AddDefaultFilterFrequencies(sweepFilter, kernelRiseTau, kernelDecayTau, PSX_FILTER_SWEEP) + PSX_AddDefaultFilterFrequencies(deconvFilter, kernelRiseTau, kernelDecayTau, PSX_FILTER_DECONV) + + parameterJsonID = JWN_GetWaveNoteAsJSON(psxKernelDataset) + parameterPath = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX + JSON_AddTreeObject(parameterJsonID, parameterPath) + JSON_AddString(parameterJsonID, parameterPath + "/id", id) + JSON_AddVariable(parameterJsonID, parameterPath + "/numberOfSDs", numberOfSDs) + JSON_AddVariable(parameterJsonID, parameterPath + "/maxTauFactor", maxTauFactor) + parameterPath = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX_RISETIME + JSON_AddTreeObject(parameterJsonID, parameterPath) + JSON_AddVariable(parameterJsonID, parameterPath + "/upperThreshold", riseTime[%$"Upper Threshold"]) + JSON_AddVariable(parameterJsonID, parameterPath + "/lowerThreshold", riseTime[%$"Lower Threshold"]) + JSON_AddVariable(parameterJsonID, parameterPath + "/differentiateThreshold", riseTime[%$"Differentiate Threshold"]) + + parameterPath = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX_DECONV_BP_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"]) + JSON_AddVariable(parameterJsonID, parameterPath + "/filterLowDefault", deconvFilter[%$"Filter Low (Default)"]) + JSON_AddVariable(parameterJsonID, parameterPath + "/filterHighDefault", deconvFilter[%$"Filter High (Default)"]) + JSON_AddVariable(parameterJsonID, parameterPath + "/filterOrderDefault", deconvFilter[%$"Filter Order (Default)"]) + + parameterPath = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX_SWEEP_BP_FILTER + JSON_AddTreeObject(parameterJsonID, parameterPath) + JSON_AddVariable(parameterJsonID, parameterPath + "/filterLow", sweepFilter[%$"Filter Low"]) + JSON_AddVariable(parameterJsonID, parameterPath + "/filterHigh", sweepFilter[%$"Filter High"]) + JSON_AddVariable(parameterJsonID, parameterPath + "/filterOrder", sweepFilter[%$"Filter Order"]) + JSON_AddVariable(parameterJsonID, parameterPath + "/filterLowDefault", sweepFilter[%$"Filter Low (Default)"]) + JSON_AddVariable(parameterJsonID, parameterPath + "/filterHighDefault", sweepFilter[%$"Filter High (Default)"]) + JSON_AddVariable(parameterJsonID, parameterPath + "/filterOrderDefault", sweepFilter[%$"Filter Order (Default)"]) + + numCombos = DimSize(psxKernelDataset, ROWS) / PSX_KERNEL_OUTPUTWAVES_PER_ENTRY + ASSERT(IsInteger(numCombos) && numCombos > 0, "Invalid number of input sets from psxKernel()") + + WAVE/WAVE output = SFH_CreateSFRefWave(exd.graph, SF_OP_PSX, numCombos * PSX_OPERATION_OUTPUT_WAVES_PER_ENTRY) + + WAVE/T labelsTemplate = ListToTextWave(PSX_EVENT_DIMENSION_LABELS, ";") + ASSERT(DimSize(labelsTemplate, ROWS) == PSX_OPERATION_OUTPUT_WAVES_PER_ENTRY, "Mismatched label wave") + Duplicate/FREE/T labelsTemplate, labels + + PSX_OperationSetDimensionLabels(output, numCombos, labels, labelsTemplate) + try - numberOfSDs = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX, 2, defValue = PSX_NUMBER_OF_SDS_DEFAULT, checkFunc = IsStrictlyPositiveAndFinite) - sweepFilterLow = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX, 3, defValue = PSX_DEFAULT_FILTER_LOW, checkFunc = IsNullOrPositiveAndFinite) - sweepFilterHigh = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX, 4, defValue = PSX_DEFAULT_FILTER_HIGH, checkFunc = IsNullOrPositiveAndFinite) - maxTauFactor = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX, 5, defValue = PSX_DEFAULT_MAX_TAU_FACTOR, checkFunc = IsStrictlyPositiveAndFinite) - WAVE riseTime = SFH_GetArgumentAsWave(exd, SF_OP_PSX, 6, defOp = "psxRiseTime()", expectedMinorType = IGOR_TYPE_64BIT_FLOAT, singleResult = 1) - ASSERT(IsNumericWave(riseTime), "Invalid return from psxRiseTime") - WAVE deconvFilter = SFH_GetArgumentAsWave(exd, SF_OP_PSX, 7, defOp = "psxDeconvFilter()", singleResult = 1) - - parameterJsonID = JWN_GetWaveNoteAsJSON(psxKernelDataset) - parameterPath = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX - JSON_AddTreeObject(parameterJsonID, parameterPath) - JSON_AddString(parameterJsonID, parameterPath + "/id", id) - JSON_AddVariable(parameterJsonID, parameterPath + "/numberOfSDs", numberOfSDs) - JSON_AddVariable(parameterJsonID, parameterPath + "/sweepFilterLow", sweepFilterLow) - JSON_AddVariable(parameterJsonID, parameterPath + "/sweepFilterHigh", sweepFilterHigh) - JSON_AddVariable(parameterJsonID, parameterPath + "/maxTauFactor", maxTauFactor) - parameterPath = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX_RISETIME - JSON_AddTreeObject(parameterJsonID, parameterPath) - JSON_AddVariable(parameterJsonID, parameterPath + "/upperThreshold", riseTime[%$"Upper Threshold"]) - JSON_AddVariable(parameterJsonID, parameterPath + "/lowerThreshold", riseTime[%$"Lower Threshold"]) - JSON_AddVariable(parameterJsonID, parameterPath + "/differentiateThreshold", riseTime[%$"Differentiate 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()") - - WAVE/WAVE output = SFH_CreateSFRefWave(exd.graph, SF_OP_PSX, numCombos * PSX_OPERATION_OUTPUT_WAVES_PER_ENTRY) + for(;;) + parameterPath = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX_DECONV_BP_FILTER + JSON_SetVariable(parameterJsonID, parameterPath + "/filterOrder", deconvFilter[%$"Filter Order"]) - path = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX_KERNEL - kernelAmp = JWN_GetNumberFromWaveNote(psxKernelDataset, path + "/amp") - ASSERT(IsFinite(kernelAmp), "psxKernel amplitude must be finite") - kernelRiseTau = JWN_GetNumberFromWaveNote(psxKernelDataset, path + "/riseTau") - ASSERT(IsFinite(kernelRiseTau), "riseTau must be finite") - kernelDecayTau = JWN_GetNumberFromWaveNote(psxKernelDataset, path + "/decayTau") - ASSERT(IsFinite(kernelDecayTau), "decayTau must be finite") + parameterPath = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX_SWEEP_BP_FILTER + JSON_SetVariable(parameterJsonID, parameterPath + "/filterOrder", sweepFilter[%$"Filter Order"]) - WAVE/T labelsTemplate = ListToTextWave(PSX_EVENT_DIMENSION_LABELS, ";") - ASSERT(DimSize(labelsTemplate, ROWS) == PSX_OPERATION_OUTPUT_WAVES_PER_ENTRY, "Mismatched label wave") - Duplicate/FREE/T labelsTemplate, labels + sprintf msg, "sweepFilter: %s, deconvFilter: %s", NumericWaveToList(sweepFilter, ";", trailSep = 0), NumericWaveToList(deconvFilter, ";", trailSep = 0) + DEBUGPRINT(msg) - PSX_OperationSetDimensionLabels(output, numCombos, labels, labelsTemplate) + writeIndex = 0 + for(i = 0; i < numCombos; i += 1) + readIndex = i + success = !PSX_OperationSweepGathering(exd.graph, psxKernelDataset, parameterJsonID, sweepFilter, deconvFilter, readIndex, writeIndex, output) + writeIndex += success + endfor - for(i = 0; i < numCombos; i += 1) - success = !PSX_OperationSweepGathering(exd.graph, psxKernelDataset, parameterJsonID, sweepFilterLow, sweepFilterHigh, deconvFilter, idx, output) - idx += success - endfor + numCombos = writeIndex - numCombos = idx + if(numCombos == 0) + Abort + endif - if(numCombos == 0) - Abort - endif + [constFilterOrderSweep, minFilterOrderSweep] = PSX_HasConstantFilterOrder(output, numCombos, PSX_FILTER_SWEEP) + sprintf msg, "constFilterOrderSweep=%g, minFilterOrderSweep=%g", constFilterOrderSweep, minFilterOrderSweep + DEBUGPRINT(msg) - Redimension/N=(numCombos * PSX_OPERATION_OUTPUT_WAVES_PER_ENTRY) output + [constFilterOrderDeconv, minFilterOrderDeconv] = PSX_HasConstantFilterOrder(output, numCombos, PSX_FILTER_DECONV) + sprintf msg, "constFilterOrderDeconv=%g, minFilterOrderDeconv=%g", constFilterOrderDeconv, minFilterOrderDeconv + DEBUGPRINT(msg) + + if(constFilterOrderSweep && constFilterOrderDeconv) + Redimension/N=(numCombos * PSX_OPERATION_OUTPUT_WAVES_PER_ENTRY) output + break + endif + + if(!constFilterOrderSweep) + sweepFilter[%$"Filter Order"] = minFilterOrderSweep + endif + + if(!constFilterOrderDeconv) + deconvFilter[%$"Filter Order"] = minFilterOrderDeconv + endif + endfor [WAVE hist, WAVE fit, peakThresh, dataUnit] = PSX_CalculatePeakThreshold(output, numCombos, numberOfSDs) WaveClear hist, fit + path = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX_KERNEL + kernelAmp = JWN_GetNumberFromWaveNote(psxKernelDataset, path + "/amp") + ASSERT(IsFinite(kernelAmp), "psxKernel amplitude must be finite") + for(i = 0; i < numCombos; i += 1) PSX_OperationImpl(exd.graph, parameterJsonID, id, peakThresh, maxTauFactor, riseTime, kernelAmp, kernelRiseTau, kernelDecayTau, i, output) endfor @@ -4732,6 +5336,9 @@ Function/WAVE PSX_Operation(STRUCT SF_ExecutionData &exd) SFH_FATAL_ERROR("Could not gather sweep data for psx") endtry + parameterPath = SF_META_USER_GROUP + PSX_JWN_PARAMETERS + "/" + SF_OP_PSX + JSON_AddVariable(parameterJsonID, parameterPath + "/peakThreshold", peakThresh) + JWN_SetWaveNoteFromJSON(output, parameterJsonID) JWN_SetStringInWaveNote(output, SF_META_DATATYPE, SF_DATATYPE_PSX) JWN_SetStringInWaveNote(output, SF_META_OPSTACK, AddListItem(SF_OP_PSX, "")) @@ -4741,6 +5348,40 @@ Function/WAVE PSX_Operation(STRUCT SF_ExecutionData &exd) return SFH_GetOutputForExecutor(output, exd.graph, SF_OP_PSX) End +Function [variable constFilterOrder, variable minFilterOrder] PSX_HasConstantFilterOrder(WAVE/WAVE output, variable numCombos, variable type) + + variable i, realOrder + string lbl, key, msg + + switch(type) + case PSX_FILTER_SWEEP: + lbl = "sweepDataOffFilt" + break + case PSX_FILTER_DECONV: + lbl = "sweepDataOffFiltDeconv" + break + default: + FATAL_ERROR("Unexpected type") + endswitch + + Make/FREE/N=(numCombos) filterOrder + + for(i = 0; i < numCombos; i += 1) + key = PSX_GenerateKey(lbl, i) + WAVE/Z wv = output[%$key] + ASSERT(WaveExists(wv), "Expected data wave") + + realOrder = JWN_GetNumberFromWaveNote(wv, "realOrder") + ASSERT(IsInteger(realOrder), "Expected an integer") + filterOrder[i] = realOrder + endfor + + sprintf msg, "filterOrder: %s", NumericWaveToList(filterOrder, ";", trailSep = 0) + DEBUGPRINT(msg) + + return [IsConstant(filterOrder, filterOrder[0]), WaveMin(filterOrder)] +End + /// @brief Implementation of the `psxKernel` operation /// // Returns a SweepFormula dataset with n * PSX_KERNEL_OUTPUTWAVES_PER_ENTRY @@ -4755,8 +5396,8 @@ End // psxKernel([select(...), riseTau, decayTau, amp]) Function/WAVE PSX_OperationKernel(STRUCT SF_ExecutionData &exd) - variable riseTau, decayTau, amp, dt, numPoints, numCombos, i, offset, idx - string parameterPath, key + variable riseTau, decayTau, amp, dt, numPoints, numCombos, i, offset, idx, shorterPerc + string parameterPath, key, bsPanel, comboKey, msg SFH_CheckArgumentCount(exd, SF_OP_PSX_KERNEL, 0, maxArgs = 4) @@ -4782,11 +5423,13 @@ Function/WAVE PSX_OperationKernel(STRUCT SF_ExecutionData &exd) Make/FREE/N=(0) allResolvedRanges Make/FREE/N=(0)/T allSelectHashes + bsPanel = BSP_GetPanel(exd.graph) + for(i = 0; i < numCombos; i += 1) - WAVE/Z sweepData = sweepDataRef[i] - ASSERT(WaveExists(sweepData), "Can't handle invalid sweepData waves") + WAVE/Z sweepDataRaw = sweepDataRef[i] + ASSERT(WaveExists(sweepDataRaw), "Can't handle invalid sweepData waves") - [WAVE singleSelectData, WAVE range] = SFH_ParseToSelectDataWaveAndRange(sweepData) + [WAVE singleSelectData, WAVE range] = SFH_ParseToSelectDataWaveAndRange(sweepDataRaw) [WAVE resolvedRanges, WAVE/T epochRangeNames] = SFH_GetNumericRangeFromEpochFromSingleSelect(exd.graph, singleSelectData, range) @@ -4794,12 +5437,24 @@ Function/WAVE PSX_OperationKernel(STRUCT SF_ExecutionData &exd) continue endif - numPoints = DimSize(sweepData, ROWS) - dt = DimDelta(sweepData, ROWS) + numPoints = DimSize(sweepDataRaw, ROWS) if(IsOdd(numPoints)) // throw away one point so that FFT works - Redimension/N=(--numPoints) sweepData + Redimension/N=(--numPoints) sweepDataRaw + endif + + WAVE sweepData = ShortenWaveForFFTIfRequired(sweepDataRaw) + + numPoints = DimSize(sweepData, ROWS) + dt = DimDelta(sweepData, ROWS) + + if(numPoints != DimSize(sweepDataRaw, ROWS)) + comboKey = PSX_GenerateComboKey(exd.graph, singleSelectData, range) + shorterPerc = (DimSize(sweepDataRaw, ROWS) - numPoints) / DimSize(sweepDataRaw, ROWS) * ONE_TO_PERCENT + sprintf msg, "psxKernel: The sweep for combination \"%s\" has a very unfitting length and was therefore shortened by %.3g%%.", comboKey, shorterPerc + SF_SetOutputState(msg, SF_MSG_WARN) + SF_DisplayOutputStateInGUI(bsPanel) endif WAVE/WAVE result = PSX_GetPSXKernel(riseTau, decayTau, amp, numPoints, dt, range) @@ -4865,38 +5520,82 @@ Function/WAVE PSX_OperationRiseTime(STRUCT SF_ExecutionData &exd) return SFH_GetOutputForExecutor(output, exd.graph, SF_OP_PSX_RISETIME) End -// psxDeconvFilter([low, high, order]) -Function/WAVE PSX_OperationDeconvFilter(STRUCT SF_ExecutionData &exd) +// psxDeconvBPFilter([low, high, order]) +Function/WAVE PSX_OperationDeconvBPFilter(STRUCT SF_ExecutionData &exd) - variable low, high, order + variable low, high, first, second, order - SFH_CheckArgumentCount(exd, SF_OP_PSX_DECONV_FILTER, 0, maxArgs = 3) + SFH_CheckArgumentCount(exd, SF_OP_PSX_DECONV_BP_FILTER, 0, maxArgs = 3) - low = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX_DECONV_FILTER, 0, defValue = NaN, checkFunc = IsNullOrPositiveAndFinite, checkDefault = 0) - high = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX_DECONV_FILTER, 1, defValue = NaN, checkFunc = IsNullOrPositiveAndFinite, checkDefault = 0) - order = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX_DECONV_FILTER, 2, defValue = NaN, checkFunc = IsStrictlyPositiveAndFinite, checkDefault = 0) + first = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX_DECONV_BP_FILTER, 0, defValue = NaN, checkFunc = IsNullOrPositiveAndFinite, checkDefault = 0) + second = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX_DECONV_BP_FILTER, 1, defValue = NaN, checkFunc = IsNullOrPositiveAndFinite, checkDefault = 0) + order = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX_DECONV_BP_FILTER, 2, defValue = NaN, checkFunc = IsStrictlyPositiveAndFinite, checkDefault = 0) - Make/D/FREE params = {low, high, order} - SetDimensionLabels(params, "Filter Low;Filter High;Filter Order", ROWS) + if(!IsNaN(first) && !IsNaN(second)) + [high, low] = MinMax(first, second) + else + low = first + high = second + endif + + if(IsOdd(order)) + order += 1 + endif - WAVE/WAVE output = SFH_CreateSFRefWave(exd.graph, SF_OP_PSX_DECONV_FILTER, 1) + Make/D/FREE params = {low, high, order, NaN, NaN, PSX_DECONV_FILTER_DEF_ORDER} + SetDimensionLabels(params, "Filter Low;Filter High;Filter Order;Filter Low (Default);Filter High (Default);Filter Order (Default)", ROWS) + + WAVE/WAVE output = SFH_CreateSFRefWave(exd.graph, SF_OP_PSX_DECONV_BP_FILTER, 1) output[0] = params - return SFH_GetOutputForExecutor(output, exd.graph, SF_OP_PSX_DECONV_FILTER) + return SFH_GetOutputForExecutor(output, exd.graph, SF_OP_PSX_DECONV_BP_FILTER) +End + +// psxSweepBPFilter([low, high, order]) +Function/WAVE PSX_OperationSweepBPFilter(STRUCT SF_ExecutionData &exd) + + variable low, high, first, second, order + + SFH_CheckArgumentCount(exd, SF_OP_PSX_SWEEP_BP_FILTER, 0, maxArgs = 3) + + first = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX_SWEEP_BP_FILTER, 0, defValue = NaN, checkFunc = IsNullOrPositiveAndFinite, checkDefault = 0) + second = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX_SWEEP_BP_FILTER, 1, defValue = NaN, checkFunc = IsNullOrPositiveAndFinite, checkDefault = 0) + order = SFH_GetArgumentAsNumeric(exd, SF_OP_PSX_SWEEP_BP_FILTER, 2, defValue = NaN, checkFunc = IsStrictlyPositiveAndFinite, checkDefault = 0) + + if(!IsNaN(first) && !IsNaN(second)) + [high, low] = MinMax(first, second) + else + low = first + high = second + endif + + if(IsOdd(order)) + order += 1 + endif + + Make/D/FREE params = {low, high, order, NaN, NaN, PSX_SWEEP_FILTER_DEF_ORDER} + SetDimensionLabels(params, "Filter Low;Filter High;Filter Order;Filter Low (Default);Filter High (Default);Filter Order (Default)", ROWS) + + WAVE/WAVE output = SFH_CreateSFRefWave(exd.graph, SF_OP_PSX_SWEEP_BP_FILTER, 1) + + output[0] = params + + return SFH_GetOutputForExecutor(output, exd.graph, SF_OP_PSX_SWEEP_BP_FILTER) End static Function/WAVE PSX_GetAllStatsProperties() - Make/FREE/T allProps = {"amp", \ - "peak", "peaktime", \ - "deconvpeak", "deconvpeaktime", \ - "baseline", "baselinetime", \ - "xinterval", \ - "tau", \ - "estate", "fstate", "fitresult", \ - "slewrate", "slewratetime", \ - "risetime", "onsettime"} + Make/FREE/T allProps = {"amp", \ + "peak", "peaktime", \ + "deconvpeak", "deconvpeaktime", \ + "baseline", "baselinetime", \ + "xinterval", \ + "slowtau", "fasttau", "weightedtau", \ + "estate", "fstate", "fitresult", \ + "slewrate", "slewratetime", \ + "onsettime", "onset", \ + "risetime", "rise"} return allProps End @@ -5325,10 +6024,7 @@ static Function PSX_SetCombo(string win, variable comboIndex) ReplaceWave/W=$extSingleGraph allinCDF SetDataFolder currentDFR - WAVE eventLocationTicks = GetPSXEventLocationTicks(comboDFR) - WAVE eventLocationLabels = GetPSXEventLocationLabels(comboDFR) - - PSX_ApplySpecialPlotProperties(psxGraph, eventLocationTicks, eventLocationLabels) + PSX_ApplySpecialPlotProperties(psxGraph) if(PSX_GetRestrictEventsToCurrentCombo(win)) PSX_UpdateAllEventGraph(win, forceAverageUpdate = 1, forceBlockIndexUpdate = 1) @@ -5545,7 +6241,7 @@ Function PSX_CheckboxProcChangeRestrictCurrentCombo(STRUCT WMCheckboxAction &cba endswitch End -Function PSX_CheckboxProcFitAcceptAverage(STRUCT WMCheckboxAction &cba) : CheckboxControl +Function PSX_CheckboxProcFitAverage(STRUCT WMCheckboxAction &cba) : CheckboxControl switch(cba.eventCode) case 2: // mouse up @@ -5556,6 +6252,19 @@ Function PSX_CheckboxProcFitAcceptAverage(STRUCT WMCheckboxAction &cba) : Checkb endswitch End +Function PSX_FitStartAmplitude(STRUCT WMSetVariableAction &sva) : SetVariableControl + + switch(sva.eventCode) + case 1: // fallthrough, mouse up + case 2: // fallthrough, Enter key + case 3: // Live update + PSX_UpdateAllEventGraph(sva.win, forceSingleEventUpdate = 1, forceAverageUpdate = 1) + break + default: + break + endswitch +End + Function PSX_PopupMenuBlockNumber(STRUCT WMPopupAction &cba) : PopupMenuControl switch(cba.eventCode) @@ -5611,7 +6320,7 @@ Function PSX_PlotStartupSettings() // propagate changes as ResizeControls queries the guides as well DoUpdate/W=$win - Make/T/FREE infoButtons = {"button_psx_info", "button_fit_results"} + Make/T/FREE infoButtons = {"button_psx_info", "button_accept_fit_results", "button_reject_fit_results", "button_undet_fit_results", "button_all_fit_results"} WAVE/T subwindows = ListToTextWave(GetAllWindows(win), ";") for(subwin : subwindows) @@ -5647,6 +6356,9 @@ Function PSX_PlotStartupSettings() endfor // default GUI values + CheckBox checkbox_show_deconv_lines, value=1, win=$win + CheckBox checkbox_show_peak_lines, value=0, win=$win + CheckBox checkbox_show_baseline_lines, value=0, win=$win CheckBox checkbox_suppress_update, value=0, win=$win ListBox listbox_select_combo, win=$win, listWave=$"", selWave=$"", helpWave=$"", selRow=0 @@ -5661,12 +6373,15 @@ Function PSX_PlotStartupSettings() CheckBox checkbox_average_events_reject, value=0, win=$specialEventPanel CheckBox checkbox_average_events_undetermined, value=0, win=$specialEventPanel CheckBox checkbox_average_events_all, value=0, win=$specialEventPanel - CheckBox checkbox_average_events_fit, value=0, win=$specialEventPanel + CheckBox checkbox_events_fit_accept, value=0, win=$specialEventPanel + CheckBox checkbox_events_fit_reject, value=0, win=$specialEventPanel + CheckBox checkbox_events_fit_undetermined, value=0, win=$specialEventPanel + CheckBox checkbox_events_fit_all, value=0, win=$specialEventPanel CheckBox checkbox_restrict_events_to_current_combination, value=0, win=$specialEventPanel - PopupMenu popupmenu_accept_fit_function, mode=1, win=$specialEventPanel SetVariable setvar_event_block_size, value=_NUM:100, win=$specialEventPanel PopupMenu popup_block, mode=1, value="", win=$specialEventPanel, userdata($PSX_UD_NUM_BLOCKS)="1" + SetVariable setvar_fit_start_amplitude, value=_NUM:20, win=$specialEventPanel StoreCurrentPanelsResizeInfo(win) @@ -5726,3 +6441,155 @@ static Function PSX_ApplyMacroToExistingPanel(string win, string mac) KillVariables/Z V_flag KillStrings/Z S_name End + +Function PSX_UpdateVisualizationHelpers(STRUCT WMCheckboxAction &cba) : CheckBoxControl + + switch(cba.eventCode) + case 2: // mouse up + PSX_ApplySpecialPlotProperties(cba.win) + break + default: + break + endswitch + + return 0 +End + +Function PSX_CalculateOnsetTimeFromAvg(WAVE AvgEvent, variable kernelAmp, variable meanOnsetTime, variable meanPeakTime) + + duplicate/FREE AvgEvent, AvgEventDiff + differentiate AvgEventDiff + smooth 500, AvgEventDiff + wavestats/Q AvgEvent + duplicate/O AvgEventDiff, root:forDisp + smooth 500, root:forDisp + + variable eventPeak, eventPeak_t, edge, backwardEdge, level, slewrate, slewrate_t + + if(kernelAmp > 0) + eventPeak = V_max + eventPeak_t = V_maxLoc + edge = FINDLEVEL_EDGE_INCREASING + backwardEdge = FINDLEVEL_EDGE_DECREASING + elseif(kernelAmp < 0) + eventPeak = V_min + eventPeak_t = V_minLoc + edge = FINDLEVEL_EDGE_DECREASING + backwardEdge = FINDLEVEL_EDGE_INCREASING + else + FATAL_ERROR("Invalid kernel amp") + endif + + variable searchEnd = -Inf // optionally: eventPeak_t - (meanPeakTime) or similar + + wavestats/Q/R=(searchEnd, eventPeak_t) AvgEventDiff + + if(kernelAmp > 0) + slewRate = v_max + slewRate_t = v_maxloc + elseif(kernelAmp < 0) + slewRate = V_min + slewRate_t = V_minLoc + endif + + level = 0.2 * slewRate + + FindLevel/R=(eventPeak_t, searchEnd)/EDGE=(edge)/Q AvgEventDiff, level + + if(!V_flag) + return V_levelX + endif + + // --- Fallback: Use second derivative to enforce slope direction --- + duplicate/FREE AvgEventDiff, d2 + differentiate d2 + smooth 200, d2 // optional smoothing for stability + + variable closestIdx = -1 + variable closestDist = Inf + variable n = DimSize(AvgEventDiff, 0) + variable xval, diffVal, d2Val + variable i + + for(i = 0; i < n; i += 1) + xval = pnt2x(AvgEventDiff, i) + + // Constrain time range + if(xval < searchEnd || xval > eventPeak_t) + continue + endif + + diffVal = AvgEventDiff[i] + d2Val = d2[i] + + // Constrain direction of slope change at the candidate level crossing + if(edge == FINDLEVEL_EDGE_INCREASING && d2Val < 0) + continue // not increasing slope + elseif(edge == FINDLEVEL_EDGE_DECREASING && d2Val > 0) + continue // not decreasing slope + endif + + // Find closest to desired level + if(abs(diffVal - level) < closestDist) + closestDist = abs(diffVal - level) + closestIdx = i + endif + endfor + + if(closestIdx >= 0) + return pnt2x(AvgEventDiff, closestIdx) + endif + + return NaN +End + +Function GetKernelDecayTau() + + string psxPath + string ListOfwin = WinList("*", ";", "WIN:64") // Only graph windows + variable i + variable num = ItemsInList(ListOfwin, ";") + + // Search all graph windows for one with psxFolder user data + for(i = 0; i < num; i += 1) + string baseWin = StringFromList(i, ListOfwin, ";") + psxPath = GetUserData(baseWin, "", "psxFolder") + if(strlen(psxPath) > 0) + break + endif + endfor + + if(strlen(psxPath) == 0) + printf "No psxFolder user data found in any graph window.\r" + return NaN + endif + + DFREF psxDFR = $psxPath + + // Search combo_* subfolders for a valid psxEvent wave + variable j = 0 + do + string folderName = GetIndexedObjName(psxPath, 4, j) // 4 = data folder + if(strlen(folderName) == 0) + break + endif + + if(stringmatch(folderName, "combo_*")) + DFREF comboDFR = psxDFR:$folderName + if(DataFolderRefStatus(comboDFR) != 0) + WAVE/Z psxEvent = comboDFR:psxEvent + if(WaveExists(psxEvent)) + variable tau = JWN_GetNumberFromWaveNote(psxEvent, "User/Parameters/psxKernel/decayTau") + if(IsFinite(tau)) + return tau + endif + endif + endif + endif + + j += 1 + while(1) + + printf "No valid psxEvent wave with decayTau found.\r" + return NaN +End diff --git a/Packages/MIES/MIES_SweepFormula_PSX_Macro.ipf b/Packages/MIES/MIES_SweepFormula_PSX_Macro.ipf index 6cc0127eb3..eede49c38c 100644 --- a/Packages/MIES/MIES_SweepFormula_PSX_Macro.ipf +++ b/Packages/MIES/MIES_SweepFormula_PSX_Macro.ipf @@ -8,7 +8,7 @@ Window PSXPanel() : Panel PauseUpdate; Silent 1 // building window... - NewPanel/K=1/W=(1251, 700, 2391, 1532) as "SweepFormula plot from " + NewPanel/K=1/W=(27, 1010, 1281, 1490) as "SweepFormula plot from " SetDrawLayer UserBack SetDrawEnv pop DrawText 47, 475, "UI" @@ -32,16 +32,15 @@ Window PSXPanel() : Panel Button button_jump_first_undet, userdata(ResizeControlsInfo)=A"!!,EN!!#=c!!#>V!!#- none -"} - Button button_psx_info, userdata="- none -" + Button button_psx_info, title="i", help={"
"}, userdata="- none -\r\n"
 	Button button_psx_info, userdata(ResizeControlsInfo)=A"!!,BI!!#CJ!!#!!#BMJ,fQL!!#](Aon\"Qzzzzzzzzzzzzzz!!#](Aon\"Qzz"
 	ListBox listbox_select_combo, userdata(ResizeControlsInfo)+=A"zzzzzzzzzzzz!!#u:Du]k- none -"}
-	Button button_fit_results, userdata="- none -"
-	GroupBox group_event, pos={14.00, 13.00}, size={69.00, 122.00}, title="Event"
+	CheckBox checkbox_events_fit_accept, pos={181.00, 33.00}, size={14.00, 14.00}, proc=PSX_CheckboxProcFitAverage
+	CheckBox checkbox_events_fit_accept, title=""
+	CheckBox checkbox_events_fit_accept, help={"Fit the accept average and store the outcome in the results wave"}
+	CheckBox checkbox_events_fit_accept, value=0
+	Button button_fit_results_accept, pos={199.00, 32.00}, size={18.00, 16.00}, proc=PSX_CopyHelpToClipboard
+	Button button_fit_results_accept, title="i", help={"
- none -
"}
+	Button button_fit_results_accept, userdata="- none -"
+	GroupBox group_event, pos={14.00, 13.00}, size={67.00, 99.00}, title="Event"
 	GroupBox group_event, help={"Toggle the display of the event traces"}
-	SetVariable setvar_event_block_size, pos={0.00, 193.00}, size={120.00, 18.00}, bodyWidth=44, proc=PSX_SetVarBlockSize
+	SetVariable setvar_event_block_size, pos={14.00, 205.00}, size={120.00, 18.00}, bodyWidth=44, proc=PSX_SetVarBlockSize
 	SetVariable setvar_event_block_size, title="Block size [%]"
 	SetVariable setvar_event_block_size, help={"Allows to restrict the all event graph to only a percentage of the events."}
 	SetVariable setvar_event_block_size, limits={0, 100, 1}, value=_NUM:100
-	PopupMenu popup_block, pos={17.00, 218.00}, size={82.00, 19.00}, bodyWidth=50, proc=PSX_PopupMenuBlockNumber
+	PopupMenu popup_block, pos={14.00, 228.00}, size={82.00, 19.00}, bodyWidth=50, proc=PSX_PopupMenuBlockNumber
 	PopupMenu popup_block, title="Block"
 	PopupMenu popup_block, help={"Select which of the event blocks to display"}
 	PopupMenu popup_block, userdata(NumberOfBlocks)="1"
 	PopupMenu popup_block, mode=1, popvalue="", value=#"\"\""
-	PopupMenu popupmenu_event_offset, pos={136.00, 168.00}, size={53.00, 19.00}, proc=PSX_PopupMenuState
+	PopupMenu popupmenu_event_offset, pos={14.00, 157.00}, size={53.00, 19.00}, proc=PSX_PopupMenuState
 	PopupMenu popupmenu_event_offset, help={"Select the time point in x direction for aligning the single event traces in the all event graph"}
 	PopupMenu popupmenu_event_offset, mode=1, popvalue="Onset", value=#"\"Onset;Peak;Slew\""
-	DefineGuide leftMenu={FL, 0.2, FR}, horizCenter={leftMenu, 0.5, FR}
+	SetVariable setvar_fit_start_amplitude, pos={14.00, 182.00}, size={144.00, 18.00}, bodyWidth=44, proc=PSX_FitStartAmplitude
+	SetVariable setvar_fit_start_amplitude, title="Fit start amplitude"
+	SetVariable setvar_fit_start_amplitude, help={"Percentage of the amplitude used to define the fit start point."}
+	SetVariable setvar_fit_start_amplitude, limits={0, 100, 1}, value=_NUM:20
+	GroupBox group_fit, pos={172.00, 14.00}, size={54.00, 98.00}, title="Fit"
+	Button button_fit_results_reject, pos={199.00, 52.00}, size={18.00, 16.00}, proc=PSX_CopyHelpToClipboard
+	Button button_fit_results_reject, title="i", help={"
- none -
"}
+	Button button_fit_results_reject, userdata="- none -"
+	CheckBox checkbox_events_fit_reject, pos={181.00, 53.00}, size={14.00, 14.00}, proc=PSX_CheckboxProcFitAverage
+	CheckBox checkbox_events_fit_reject, title=""
+	CheckBox checkbox_events_fit_reject, help={"Fit the reject average and store the outcome in the results wave"}
+	CheckBox checkbox_events_fit_reject, value=0
+	Button button_fit_results_undetermined, pos={199.00, 70.00}, size={18.00, 16.00}, proc=PSX_CopyHelpToClipboard
+	Button button_fit_results_undetermined, title="i", help={"
- none -
"}
+	Button button_fit_results_undetermined, userdata="- none -"
+	CheckBox checkbox_events_fit_undetermined, pos={181.00, 72.00}, size={14.00, 14.00}, proc=PSX_CheckboxProcFitAverage
+	CheckBox checkbox_events_fit_undetermined, title=""
+	CheckBox checkbox_events_fit_undetermined, help={"Fit the undet average and store the outcome in the results wave"}
+	CheckBox checkbox_events_fit_undetermined, value=0
+	Button button_fit_results_all, pos={199.00, 90.00}, size={18.00, 16.00}, proc=PSX_CopyHelpToClipboard
+	Button button_fit_results_all, title="i", help={"
- none -
"}
+	Button button_fit_results_all, userdata="- none -"
+	CheckBox checkbox_events_fit_all, pos={181.00, 91.00}, size={14.00, 14.00}, proc=PSX_CheckboxProcFitAverage
+	CheckBox checkbox_events_fit_all, title=""
+	CheckBox checkbox_events_fit_all, help={"Fit the all average and store the outcome in the results wave"}
+	CheckBox checkbox_events_fit_all, value=0
+	DefineGuide leftMenu={FL, 0.214995, FR}, horizCenter={leftMenu, 0.5, FR}
 	SetWindow kwTopWin, hook(resetScaling)=IH_ResetScaling
 	SetWindow kwTopWin, hook(ctrl)=PSX_AllEventGraphHook
-	Execute/Q/Z "SetWindow kwTopWin sizeLimit={750,200.25,inf,inf}" // sizeLimit requires Igor 7 or later
+	Execute/Q/Z "SetWindow kwTopWin sizeLimit={820,200.25,inf,inf}" // sizeLimit requires Igor 7 or later
 	Display/W=(225, 62, 675, 188)/FG=(horizCenter, FT, FR, FB)/HOST=#
 	RenameWindow #, Single
 	SetActiveSubwindow ##
-	Display/W=(225, 62, 675, 188)/FG=(leftMenu, FT, horizCenter, FB)/HOST=#
+	Display/W=(287, 62, 675, 188)/FG=(leftMenu, FT, horizCenter, FB)/HOST=#
 	RenameWindow #, All
 	SetActiveSubwindow ##
 	RenameWindow #, SpecialEventPanel
diff --git a/Packages/MIES/MIES_Utilities_Algorithm.ipf b/Packages/MIES/MIES_Utilities_Algorithm.ipf
index 9f9c49379f..64318447a5 100644
--- a/Packages/MIES/MIES_Utilities_Algorithm.ipf
+++ b/Packages/MIES/MIES_Utilities_Algorithm.ipf
@@ -6,6 +6,8 @@
 #pragma ModuleName = MIES_UTILS_ALGORITHM
 #endif // AUTOMATED_TESTING
 
+static Constant FFT_PRIME_FACTOR_LIMIT = 30
+
 /// @file MIES_Utilities_Algorithm.ipf
 /// @brief utility functions for common algorithms
 
@@ -1040,6 +1042,36 @@ threadsafe Function DoPowerSpectrum(WAVE input, WAVE output, variable col)
 	output[][col] = magsqr(powerSpectrum[p])
 End
 
+/// @brief Cut off some trailing points from the wave to allow fast FFT
+Function/WAVE ShortenWaveForFFTIfRequired(WAVE input)
+
+	variable numRowsOpt, numRows
+	string key
+
+	numRows = DimSize(input, ROWS)
+	WAVE primes = GetPrimeFactors(numRows)
+
+	if(WaveMax(primes) > 1000)
+
+		key = CA_GetGoodFFTSizesKeys()
+		WAVE/Z goodFFTSizes = CA_TryFetchingEntryFromCache(key)
+
+		if(!WaveExists(goodFFTSizes))
+			WAVE goodFFTSizes = GetGoodFFTSizes()
+			CA_StoreEntryIntoCache(key, goodFFTSizes)
+		endif
+
+		FindLevel/Q goodFFTSizes, numRows
+		ASSERT_TS(!V_flag, "Could not find good FFT size")
+
+		numRowsOpt = goodFFTSizes[trunc(V_LevelX)]
+		Duplicate/FREE/RMD=[0, numRowsOpt - 1] input, inputOpt
+		return inputOpt
+	endif
+
+	return input
+End
+
 /// @brief Perform FFT on input with optionally given window function
 ///
 /// @param input   Wave to perform FFT on
@@ -1275,3 +1307,65 @@ Function FindSequenceReverseWrapper(WAVE sequence, WAVE source)
 		start      = foundIndex + 1
 	endfor
 End
+
+/// @brief Return all prime factors from num as free wave
+threadsafe Function/WAVE GetPrimeFactors(variable num)
+
+	PrimeFactors/Q num
+
+	WAVE W_primeFactors
+	MakeWaveFree(W_primeFactors)
+
+	return W_primeFactors
+End
+
+threadsafe static Function/WAVE FFTHelperSeriesImpl(variable n, variable a, variable b)
+
+	ASSERT_TS(n >= 1 && n < 30, "Invalid range")
+
+	variable numEntries = ceil(n / 3)
+	variable lambda, i
+
+	Make/FREE/D/N=(numEntries) results
+
+	for(i = 0; i < numEntries; i += 1)
+		lambda = n - (a + 3 * i)
+		if(lambda < 0)
+			break
+		endif
+
+		results[i] = 3^(b + 2 * i) * 2^(lambda)
+	endfor
+
+	Redimension/N=(i) results
+
+	return results
+End
+
+threadsafe Function/WAVE FFTHelperSilverSeries(variable n)
+
+	return FFTHelperSeriesImpl(n, 3, 2)
+End
+
+threadsafe Function/WAVE FFTHelperGoldSeries(variable n)
+
+	return FFTHelperSeriesImpl(n, 1, 1)
+End
+
+// @brief Returns a wave with all values which only have 2's and 3's as prime factors up to 2^30
+threadsafe Function/WAVE GetGoodFFTSizes()
+
+	Make/FREE/WAVE/N=(FFT_PRIME_FACTOR_LIMIT) gold, silver
+	gold[0]      = NewFreeWave(IGOR_TYPE_64BIT_FLOAT, 0)
+	gold[1, Inf] = FFTHelperGoldSeries(p)
+
+	silver[0]      = NewFreeWave(IGOR_TYPE_64BIT_FLOAT, 0)
+	silver[1, Inf] = FFTHelperSilverSeries(p)
+
+	Concatenate/FREE/NP=(ROWS) {gold, silver}, resultsWave
+	Concatenate/FREE/NP=(ROWS) {resultsWave}, results
+
+	Sort results, results
+
+	return results
+End
diff --git a/Packages/MIES/MIES_WaveDataFolderGetters.ipf b/Packages/MIES/MIES_WaveDataFolderGetters.ipf
index 3ae72cf33d..cd0d82b478 100644
--- a/Packages/MIES/MIES_WaveDataFolderGetters.ipf
+++ b/Packages/MIES/MIES_WaveDataFolderGetters.ipf
@@ -8387,15 +8387,15 @@ End
 /// @name SweepFormula PSX
 ///@{
 
-static Constant PSX_WAVE_VERSION       = 3
-static Constant PSX_EVENT_WAVE_COLUMNS = 17
+static Constant PSX_WAVE_VERSION       = 5
+static Constant PSX_EVENT_WAVE_COLUMNS = 21
 
 /// @brief Return the upgraded psxEvent wave
 Function/WAVE UpgradePSXEventWave(WAVE psxEvent)
 
 	if(WaveVersionIsAtLeast(psxEvent, PSX_WAVE_VERSION))
 		return psxEvent
-	elseif(WaveVersionIsAtLeast(psxEvent, 2))
+	elseif(WaveVersionIsAtLeast(psxEvent, 2)) // Version 2-4
 
 		if(!AlreadyCalledOnce(CO_PSX_UPGRADE_EVENT))
 			print "The algorithm for psp/psc event detection was heavily overhauled, therefore we are very sorry " \
@@ -8423,25 +8423,30 @@ End
 /// -  1/deconvPeak: event amplitude in deconvoluted data [y unit of data]
 /// -  2/deconvPeak_t: deconvolved peak time [ms]
 /// -  3/peak: Maximum (positive kernel amp sign) or minimum (negative kernel amp sign) in the range of
-///            [deconvPeak_t – kernelRiseTau or devonvPeak_t of the previous event (whichever comes later),
-///             deconvPeak_t + 0.33 * kernelDecayTau or deconvPeak_t of the next event (which ever comes first)]
+///            [baseline_t – 5 * kernelRiseTau or baseline_t or devonvPeak_t/peak_t of the previous event (whichever comes later),
+///             baseline_t + 0.33 * kernelDecayTau or deconvPeak_t of the next event (which ever comes first)]
 ///            in the filtered sweep wave
 /// -  4/peak_t: peak time
 /// -  5/baseline: Maximum (negative kernel amp sign) or minimum (positive kernel amp sign) in the range of
-///                    [peak_t – 10 * kernelRiseTau, peak_t], averaged over +/- 5 points, in the filtered sweep wave
+///                    [deconvPeak_t – 10 * kernelRiseTau or the minimum of peak_t and deconvPeak_t of previous event (whichever comes later),
+///                     deconvPeak_t], averaged over +/- 5 points, in the filtered sweep wave
 /// -  6/baseline_t: baseline time
 /// -  7/amplitude: Relative amplitude: [3] - [5]
 /// -  8/iei: Time difference to previous event (inter event interval) [ms]
-/// -  9/tau: Decay constant tau of exponential fit
-/// - 10/Fit manual QC call: One of @ref PSXStates
-/// - 11/Fit result: 1 for success, everything smaller than 0 is failure:
+/// -  9/weightedTau: Weighted tau of the double exponential fit
+/// - 10/fastTau: Fast tau of the double exponential fit
+/// - 11/slowTau: Slow tau of the double exponential fit
+/// - 12/Fit manual QC call: One of @ref PSXStates
+/// - 13/Fit result: 1 for success, everything smaller than 0 is failure:
 ///   - `]-10000, 0[`: CurveFit error codes
 ///   - `]-inf, -10000]`: Custom error codes, one of @ref FitEventDecayCustomErrors
-/// - 12/Event manual QC call: One of @ref PSXStates
-/// - 13/Onset time as calculated by PSX_CalculateOnsetTime
-/// - 14/Rise Time as calculated by PSX_CalculateRiseTime
-/// - 15/Slew Rate
-/// - 16/Slew Rate Time
+/// - 14/Event manual QC call: One of @ref PSXStates
+/// - 15/Onset Time as calculated by PSX_CalculateOnsetTime
+/// - 16/Onset y-value
+/// - 17/Rise Time as calculated by PSX_CalculateRiseTime
+/// - 18/Rise y-value
+/// - 19/Slew Rate
+/// - 20/Slew Rate Time
 Function/WAVE GetPSXEventWaveAsFree()
 
 	variable versionOfWave = PSX_WAVE_VERSION
@@ -8466,14 +8471,18 @@ static Function SetPSXEventDimensionLabels(WAVE wv)
 	SetDimLabel COLS, 6, baseline_t, wv
 	SetDimLabel COLS, 7, amplitude, wv
 	SetDimLabel COLS, 8, iei, wv
-	SetDimLabel COLS, 9, tau, wv
-	SetDimLabel COLS, 10, $"Fit manual QC call", wv
-	SetDimLabel COLS, 11, $"Fit result", wv
-	SetDimLabel COLS, 12, $"Event manual QC call", wv
-	SetDimLabel COLS, 13, $"Onset Time", wv
-	SetDimLabel COLS, 14, $"Rise Time", wv
-	SetDimLabel COLS, 15, $"Slew Rate", wv
-	SetDimLabel COLS, 16, $"Slew Rate Time", wv
+	SetDimLabel COLS, 9, weightedTau, wv
+	SetDimLabel COLS, 10, slowTau, wv
+	SetDimLabel COLS, 11, fastTau, wv
+	SetDimLabel COLS, 12, $"Fit manual QC call", wv
+	SetDimLabel COLS, 13, $"Fit result", wv
+	SetDimLabel COLS, 14, $"Event manual QC call", wv
+	SetDimLabel COLS, 15, $"Onset Time", wv
+	SetDimLabel COLS, 16, $"Onset", wv
+	SetDimLabel COLS, 17, $"Rise Time", wv
+	SetDimLabel COLS, 18, $"Rise", wv
+	SetDimLabel COLS, 19, $"Slew Rate", wv
+	SetDimLabel COLS, 20, $"Slew Rate Time", wv
 End
 
 Function/WAVE GetPSXSingleEventFitWaveFromDFR(DFREF dfr)
@@ -8575,14 +8584,41 @@ Function/WAVE GetPSXSweepDataOffFiltDeconvWaveFromDFR(DFREF dfr)
 	return GetWaveFromFolder(dfr, "sweepDataOffFiltDeconv")
 End
 
-Function/WAVE GetPSXEventLocationLabels(DFREF dfr)
+Function/WAVE GetPSXEventLabelsAsFree(variable numEvents, string suffix)
 
-	return GetWaveFromFolder(dfr, "eventLocationLabels")
+	string name = "eventLabels_" + CleanupName(suffix, 0)
+
+	Make/FREE/T/N=(numEvents, 2) wv
+	ChangeFreeWaveName(wv, name)
+
+	SetDimLabel COLS, 1, $"Tick Type", wv
+	wv[][1] = "Major"
+
+	return wv
+End
+
+Function/WAVE GetPSXEventLabelsFromDFR(DFREF dfr, string suffix)
+
+	string name = "eventLabels_" + CleanupName(suffix, 0)
+
+	return GetWaveFromFolder(dfr, name)
+End
+
+Function/WAVE GetPSXEventTicksAsFree(variable numEvents, string suffix)
+
+	string name = "eventTicks_" + CleanupName(suffix, 0)
+
+	Make/D/FREE/N=(numEvents) wv
+	ChangeFreeWaveName(wv, name)
+
+	return wv
 End
 
-Function/WAVE GetPSXEventLocationTicks(DFREF dfr)
+Function/WAVE GetPSXEventTicksFromDFR(DFREF dfr, string suffix)
 
-	return GetWaveFromFolder(dfr, "eventLocationTicks")
+	string name = "eventTicks_" + CleanupName(suffix, 0)
+
+	return GetWaveFromFolder(dfr, name)
 End
 
 Function/S GetPSXFolderForComboAsString(DFREF dfr, variable index)
@@ -8641,15 +8677,62 @@ Function/WAVE GetPSXAverageWave(DFREF dfr, variable state)
 	return wv
 End
 
-Function/WAVE GetPSXAcceptedAverageFitWaveFromDFR(DFREF dfr)
+Function/WAVE GetPSXRiseAverageFitWaveFromDFR(DFREF dfr, variable state)
+
+	string name = "riseAverageFit_" + PSX_StateToString(state)
+
+	WAVE/Z/D/SDFR=dfr wv = $name
+
+	if(WaveExists(wv))
+		return wv
+	endif
+
+	Make/D/N=(0) dfr:$name/WAVE=wv
+
+	return wv
+End
+
+Function/WAVE GetPSXDecayAverageFitWaveFromDFR(DFREF dfr, variable state)
 
-	WAVE/Z/D/SDFR=dfr wv = acceptedAverageFit
+	string name = "decayAverageFit_" + PSX_StateToString(state)
+
+	WAVE/Z/D/SDFR=dfr wv = $name
+
+	if(WaveExists(wv))
+		return wv
+	endif
+
+	Make/D/N=(0) dfr:$name/WAVE=wv
+
+	return wv
+End
+
+Function/WAVE GetPSXAverageFitWaveFromDFR(DFREF dfr, variable state)
+
+	string name = "averageFit_" + PSX_StateToString(state)
+
+	WAVE/Z/D/SDFR=dfr wv = $name
+
+	if(WaveExists(wv))
+		return wv
+	endif
+
+	Make/D/N=(0) dfr:$name/WAVE=wv
+
+	return wv
+End
+
+Function/WAVE GetPSXAverageFitResultsWaveFromDFR(DFREF dfr, variable state)
+
+	string name = "averageFitResults_" + PSX_StateToString(state)
+
+	WAVE/Z/T/SDFR=dfr wv = $name
 
 	if(WaveExists(wv))
 		return wv
 	endif
 
-	Make/D/N=(0) dfr:acceptedAverageFit/WAVE=wv
+	Make/T/N=(0, 2) dfr:$name/WAVE=wv
 
 	return wv
 End
diff --git a/Packages/doc/SweepFormula.rst b/Packages/doc/SweepFormula.rst
index b956b07613..e92d17cf31 100644
--- a/Packages/doc/SweepFormula.rst
+++ b/Packages/doc/SweepFormula.rst
@@ -1278,7 +1278,7 @@ The `psx` operation allows to classify miniature PSC/PSP's interactively.
 
 .. code-block:: bash
 
-   psx(id, [psxKernel(), numSDs, filterLow, filterHigh, maxTauFactor, psxRiseTime(), psxDeconvFilter()])
+   psx(id, [psxKernel(), numSDs, filterLow, filterHigh, maxTauFactor, psxRiseTime(), psxDeconvBPFilter()])
 
 The function accepts one to seven arguments.
 
@@ -1292,11 +1292,8 @@ psxKernel
 numSDs
   Number of standard deviations for the gaussian fit of the all points histogram, defaults to 2.5
 
-filterLow
-  low threshold for the bandpass filter, defaults to 550 Hz
-
-filterHigh
-  high threshold for the bandpass filter, defaults to 0 Hz
+psxSweepBPFilter
+  results from the `psxSweepBPFilter` operation
 
 maxTauFactor
   maximum tau factor, the decay tau from fitting the event must be smaller than the fit range
@@ -1305,12 +1302,22 @@ maxTauFactor
 psxRiseTime
   results from the `psxRiseTime` operation
 
-psxDeconvFilter
-  results from the `psxDeconvFilter` operation
+psxDeconvBPFilter
+  results from the `psxDeconvBPFilter` operation
 
 The plotting is implemented in a custom way. Due to that multiple `psx`
 operations can only be separated by `with` and not `and`.
 
+The filter order is internally made even as there is no difference in filter
+order `n` and `n + 1` due to implementation details of the used operation
+`FilterIIR`.
+
+The filtering for both the sweep data and the deconvoluted data uses a backing down
+algorithm for determining the filter order. The implementation starts with the
+given order and decrements it by two as long as the filtering is not
+successfull for all sweeps. If we reach zero we bail out. The used filter order
+is stored in the wave note.
+
 .. code-block:: bash
 
    psx(myID)
@@ -1335,10 +1342,10 @@ select
   selections and range to operate on from the `select` operation
 
 riseTau
-  Time constant for kernel, defaults to 1
+  Time constant for kernel, defaults to 1ms
 
 decayTau
-  Time constant for kernel, defaults to 15
+  Time constant for kernel, defaults to 15ms
 
 amp
    Amplitude for kernel, defaults to -5
@@ -1393,13 +1400,43 @@ diffThreshold
    psxRiseTime(0.5, 0.9)
    psxRiseTime(0.5, 0.9, 0.15)
 
-psxDeconvFilter
-"""""""""""""""
+psxDeconvBPFilter
+"""""""""""""""""
+
+The `psxDeconvBPFilter` operation is a helper operation for `psx` to manage the deconvolution filter settings.
+This filter is a bandpass filter.
+
+   psxDeconvBPFilter([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
+`7`, for the frequencies `500` (`lowFreq`) and `50` (`highFreq`).
+Here `lowFreq` is the end and `highFreq` the start of the
+passband, see also the description of `/LO` and `/HI` from `FilterIIR`. If the
+frequency values are not ordered correctly, they are swapped.
+
+.. code-block:: bash
+
+   psxDeconvBPFilter(800, 100)
+   psxDeconvBPFilter(400, 50, 11)
+
+psxSweepBPFilter
+"""""""""""""""""
 
-The `psxDeconvFilter` operation is a helper operation for `psx` to manage the deconvolution filter settings.
+The `psxSweepBPFilter` operation is a helper operation for `psx` to manage the sweep filter settings.
 This filter is a bandpass filter.
 
-   psxDeconvFilter([lowFreq, highFreq, order])
+   psxSweepBPFilter([lowFreq, highFreq, order])
 
 The function accepts zero to three arguments.
 
@@ -1415,12 +1452,13 @@ order
 The default values of `NaN` are replaced inside `psx`. For the order this is
 `7`, for the frequencies `500` (`lowFreq`) and `50` (`highFreq`).
 Here `lowFreq` is the end and `highFreq` the start of the
-passband, see also the description of `/LO` and `/HI` from `FilterIIR`.
+passband, see also the description of `/LO` and `/HI` from `FilterIIR`. If the
+frequency values are not ordered correctly, they are swapped.
 
 .. code-block:: bash
 
-   psxDeconvFilter(800, 100)
-   psxDeconvFilter(400, 50, 11)
+   psxDeconvBPFilter(800, 100)
+   psxDeconvBPFilter(400, 50, 11)
 
 psxstats
 """"""""
@@ -1453,7 +1491,7 @@ select
 prop
   column of the `psx` event results waves to plot.
   Choices are: `amp`, `peak`, `peaktime`, `deconvpeak`, `deconvpeaktime`, `baseline`, `baselinetime`, `xinterval`,
-  `tau`, `estate`, `fstate`, `fitresult`, `slewrate`, `slewratetime`, `risetime`, `onsettime`
+  `slowtau`, `fasttau`, `weightedtau`, `estate`, `fstate`, `fitresult`, `slewrate`, `slewratetime`, `risetime`, `onsettime`
 
 state
   QC state to select the events.
diff --git a/Packages/tests/Basic/UTF_SweepFormula_PSX.ipf b/Packages/tests/Basic/UTF_SweepFormula_PSX.ipf
index 03846f0cf9..2ea18602c8 100644
--- a/Packages/tests/Basic/UTF_SweepFormula_PSX.ipf
+++ b/Packages/tests/Basic/UTF_SweepFormula_PSX.ipf
@@ -1352,9 +1352,9 @@ static Function CheckEventDataHelper(WAVE/Z/WAVE dataWref, variable index, varia
 
 	INFO("index = %d, V_numNaNs = %d, kernelAmpSign = %d", n0 = index, n1 = V_numNans, n2 = kernelAmpSign)
 
-	// 1 NaN for the first event only, the rest is onset Time
+	// 5 NaNs for the first event only, the rest is onset Time
 	if(kernelAmpSign == 1)
-		CHECK_EQUAL_VAR(V_numNaNs, 1)
+		CHECK_EQUAL_VAR(V_numNaNs, 5)
 	elseif(kernelAmpSign == -1)
 		CHECK_EQUAL_VAR(V_numNaNs, 9)
 	else
@@ -1696,6 +1696,50 @@ static Function MouseSelectionStatsPostProcNonFinite()
 	CheckPSXEventField({psxEvent_1}, {"Fit manual QC call", "Event manual QC call"}, {0}, PSX_ACCEPT)
 End
 
+static Function KeyboardSelectionPSXStatsVS()
+
+	string win, browser, code, psxGraph, psxStatsGraph, postProc, comboKey, tracenames, trace
+	variable numEvents, logMode
+
+	Make/FREE/T/N=2 combos
+	sprintf comboKey, "Range[50, 150], Sweep [0], Channel [AD6], Device [ITC16_Dev_0], Experiment [%s]", GetExperimentName()
+	combos[0] = comboKey
+	sprintf comboKey, "Range[50, 150], Sweep [2], Channel [AD6], Device [ITC16_Dev_0], Experiment [%s]", GetExperimentName()
+	combos[1] = comboKey
+	WAVE overrideResults = MIES_PSX#PSX_CreateOverrideResults(4, combos)
+
+	overrideResults[][][%$"Fit Result"]      = 1
+	overrideResults[][][%$"Tau"]             = 1
+	overrideResults[][][%$"KernelAmpSignQC"] = 1
+
+	browser = SetupDatabrowserWithSomeData()
+
+	code = GetTestCode("nothing") + "\n vs [1, 2]"
+
+	// combo0 is the current one
+
+	win = ExecuteSweepFormulaCode(browser, code)
+
+	psxGraph      = MIES_PSX#PSX_GetPSXGraph(win)
+	psxStatsGraph = MIES_PSX#PSX_GetPSXStatsGraph(psxGraph)
+
+	REQUIRE(WindowExists(psxStatsGraph))
+
+	SetActiveSubwindow $psxStatsGraph
+
+	tracenames = TraceNameList(psxStatsGraph, ";", 1)
+	trace      = StringFromList(0, tracenames)
+	CHECK_PROPER_STR(trace)
+
+	Cursor/W=$psxStatsGraph/P A, $trace, 0
+
+	CheckCurrentEvent(psxStatsGraph, 0, 0, 0)
+
+	SendKey(psxGraph, LEFT_KEY)
+	SendKey(psxGraph, LEFT_KEY)
+	CHECK_NO_RTE()
+End
+
 static Function/WAVE GetTracesHelper(string win, variable options)
 
 	return ListToTextWave(SortList(TraceNameList(win, ";", options)), ";")
@@ -2815,6 +2859,9 @@ static Function KeyboardInteractionsStatsPostProcNonFinite()
 	overrideResults[1][%$combos[0]][%$"Fit Result"] = 1
 	overrideResults[1][%$combos[0]][%$"Tau"]        = -Inf
 
+	overrideResults[3][%$combos[0]][%$"Fit Result"] = 1
+	overrideResults[3][%$combos[0]][%$"Tau"]        = -Inf
+
 	overrideResults[0][%$combos[0]][%$"Fit Result"] = 0
 	overrideResults[0][%$combos[0]][%$"Tau"]        = NaN
 
@@ -2847,11 +2894,10 @@ static Function KeyboardInteractionsStatsPostProcNonFinite()
 
 	[WAVE psxEvent_0, WAVE psxEvent_1] = GetPSXEventWavesHelper(psxStatsGraph)
 
-	CheckPSXEventField({psxEvent_0}, {"Fit manual QC call"}, {0, 3}, PSX_REJECT)
-	CheckPSXEventField({psxEvent_0}, {"Fit manual QC call"}, {1, 2}, PSX_UNDET)
+	CheckPSXEventField({psxEvent_0}, {"Fit manual QC call"}, {0}, PSX_REJECT)
+	CheckPSXEventField({psxEvent_0}, {"Fit manual QC call"}, {1, 2, 3}, PSX_UNDET)
 	CheckPSXEventField({psxEvent_0}, {"Event manual QC call"}, {0, 1, 2, 3}, PSX_UNDET)
-	CheckPSXEventField({psxEvent_1}, {"Fit manual QC call"}, {0, 1}, PSX_UNDET)
-	CheckPSXEventField({psxEvent_1}, {"Fit manual QC call"}, {2}, PSX_REJECT)
+	CheckPSXEventField({psxEvent_1}, {"Fit manual QC call"}, {0, 1, 2}, PSX_UNDET)
 	CheckPSXEventField({psxEvent_1}, {"Event manual QC call"}, {0, 1, 2}, PSX_UNDET)
 
 	SendKey(psxStatsGraph, UP_KEY)
@@ -2859,7 +2905,7 @@ static Function KeyboardInteractionsStatsPostProcNonFinite()
 	SendKey(psxStatsGraph, DOWN_KEY)
 	SendKey(psxStatsGraph, UP_KEY)
 
-	CheckCurrentEvent(psxStatsGraph, 1, 2, 4)
+	CheckCurrentEvent(psxStatsGraph, 1, 0, 3)
 
 	CheckPSXEventField({psxEvent_0}, {"Fit manual QC call"}, {0, 3}, PSX_REJECT)
 	CheckPSXEventField({psxEvent_0}, {"Fit manual QC call"}, {2}, PSX_UNDET)
@@ -2871,8 +2917,7 @@ static Function KeyboardInteractionsStatsPostProcNonFinite()
 	CheckPSXEventField({psxEvent_0}, {"Event manual QC call"}, {3}, PSX_REJECT)
 
 	CheckPSXEventField({psxEvent_1}, {"Fit manual QC call"}, {0}, PSX_ACCEPT)
-	CheckPSXEventField({psxEvent_1}, {"Fit manual QC call"}, {1}, PSX_UNDET)
-	CheckPSXEventField({psxEvent_1}, {"Fit manual QC call"}, {2}, PSX_REJECT)
+	CheckPSXEventField({psxEvent_1}, {"Fit manual QC call"}, {1, 2}, PSX_UNDET)
 	CheckPSXEventField({psxEvent_1}, {"Event manual QC call"}, {0}, PSX_ACCEPT)
 	CheckPSXEventField({psxEvent_1}, {"Event manual QC call"}, {1}, PSX_UNDET)
 	CheckPSXEventField({psxEvent_1}, {"Event manual QC call"}, {2}, PSX_UNDET)
@@ -3395,44 +3440,46 @@ static Function [variable filterLow, variable filterHigh, variable filterOrder]
 	return [filterLow, filterHigh, filterOrder]
 End
 
-static Function TestOperationDeconvFilter()
+static Function TestOperationDeconvBPFilter()
 
 	string win, str
 	variable filterLow, filterHigh, filterOrder
 
 	win = SetupDatabrowserWithSomeData()
 
-	str = "psxDeconvFilter()"
+	str = "psxDeconvBPFilter()"
 	WAVE/WAVE dataWref = SFE_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)"
+	str = "psxDeconvBPFilter(40)"
 	WAVE/WAVE dataWref = SFE_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)"
+	str = "psxDeconvBPFilter(40, 50)"
 	WAVE/WAVE dataWref = SFE_ExecuteFormula(str, win, useVariables = 0)
 	[filterLow, filterHigh, filterOrder] = TestDevonvFilterContainer(dataWref)
-	CHECK_EQUAL_VAR(filterLow, 40)
-	CHECK_EQUAL_VAR(filterHigh, 50)
+	// we automatically fixup the order for the user
+	CHECK_EQUAL_VAR(filterLow, 50)
+	CHECK_EQUAL_VAR(filterHigh, 40)
 	CHECK_EQUAL_VAR(filterOrder, NaN)
 
-	str = "psxDeconvFilter(40, 50, 11)"
+	str = "psxDeconvBPFilter(40, 50, 11)"
 	WAVE/WAVE dataWref = SFE_ExecuteFormula(str, win, useVariables = 0)
 	[filterLow, filterHigh, filterOrder] = TestDevonvFilterContainer(dataWref)
-	CHECK_EQUAL_VAR(filterLow, 40)
-	CHECK_EQUAL_VAR(filterHigh, 50)
+	// we automatically fixup the order for the user
+	CHECK_EQUAL_VAR(filterLow, 50)
+	CHECK_EQUAL_VAR(filterHigh, 40)
 	CHECK_EQUAL_VAR(filterOrder, 11)
 
 	// check parameters
 	try
-		str = "psxDeconvFilter(-1)"
+		str = "psxDeconvBPFilter(-1)"
 		WAVE/WAVE dataWref = SFE_ExecuteFormula(str, win, useVariables = 0)
 		FAIL()
 	catch
@@ -3440,7 +3487,7 @@ static Function TestOperationDeconvFilter()
 	endtry
 
 	try
-		str = "psxDeconvFilter(1, -1)"
+		str = "psxDeconvBPFilter(1, -1)"
 		WAVE/WAVE dataWref = SFE_ExecuteFormula(str, win, useVariables = 0)
 		FAIL()
 	catch
@@ -3448,7 +3495,7 @@ static Function TestOperationDeconvFilter()
 	endtry
 
 	try
-		str = "psxDeconvFilter(1, 1, -1)"
+		str = "psxDeconvBPFilter(1, 1, -1)"
 		WAVE/WAVE dataWref = SFE_ExecuteFormula(str, win, useVariables = 0)
 		FAIL()
 	catch
diff --git a/Packages/tests/Basic/UTF_Utils_Algorithm.ipf b/Packages/tests/Basic/UTF_Utils_Algorithm.ipf
index cfcf370054..237f50f50a 100644
--- a/Packages/tests/Basic/UTF_Utils_Algorithm.ipf
+++ b/Packages/tests/Basic/UTF_Utils_Algorithm.ipf
@@ -1326,3 +1326,93 @@ static Function TestFindSequenceReverseWrapper()
 End
 
 /// @}
+
+static Function TestFFTHelperSilverSeries()
+
+	try
+		MIES_UTILS_ALGORITHM#FFTHelperSilverSeries(0)
+		FAIL()
+	catch
+		CHECK_NO_RTE()
+	endtry
+
+	try
+		MIES_UTILS_ALGORITHM#FFTHelperSilverSeries(30)
+		FAIL()
+	catch
+		CHECK_NO_RTE()
+	endtry
+
+	WAVE results = MIES_UTILS_ALGORITHM#FFTHelperSilverSeries(4)
+	CHECK_EQUAL_WAVES(results, {18}, mode = WAVE_DATA)
+
+	WAVE results = MIES_UTILS_ALGORITHM#FFTHelperSilverSeries(5)
+	CHECK_EQUAL_WAVES(results, {36}, mode = WAVE_DATA)
+
+	WAVE results = MIES_UTILS_ALGORITHM#FFTHelperSilverSeries(7)
+	CHECK_EQUAL_WAVES(results, {144, 162}, mode = WAVE_DATA)
+End
+
+static Function TestFFTHelperGoldSeries()
+
+	try
+		MIES_UTILS_ALGORITHM#FFTHelperGoldSeries(0)
+		FAIL()
+	catch
+		CHECK_NO_RTE()
+	endtry
+
+	try
+		MIES_UTILS_ALGORITHM#FFTHelperGoldSeries(30)
+		FAIL()
+	catch
+		CHECK_NO_RTE()
+	endtry
+
+	WAVE results = MIES_UTILS_ALGORITHM#FFTHelperGoldSeries(4)
+	CHECK_EQUAL_WAVES(results, {24, 27}, mode = WAVE_DATA)
+
+	WAVE results = MIES_UTILS_ALGORITHM#FFTHelperGoldSeries(5)
+	CHECK_EQUAL_WAVES(results, {48, 54}, mode = WAVE_DATA)
+
+	WAVE results = MIES_UTILS_ALGORITHM#FFTHelperGoldSeries(7)
+	CHECK_EQUAL_WAVES(results, {192, 216, 243}, mode = WAVE_DATA)
+End
+
+static Function CheckPrimeFactors(variable num)
+
+	variable numEntries
+
+	WAVE   primes = GetPrimeFactors(num)
+	WAVE/Z result = GetUniqueEntries(primes)
+	CHECK_WAVE(result, NUMERIC_WAVE)
+	numEntries = DimSize(result, ROWS)
+
+	Sort result, result
+
+	INFO("result: %s", s0 = NumericWaveToList(result, ", ", trailSep = 0))
+
+	switch(numEntries)
+		case 1:
+			CHECK(result[0] == 2 || result[0] == 3) // NOLINT
+			break
+		case 2:
+			CHECK_EQUAL_WAVES(result, {2, 3}, mode = WAVE_DATA)
+			break
+		default:
+			FAIL()
+	endswitch
+End
+
+static Function TestGetGoodFFTSizes()
+
+	variable numEntries
+
+	WAVE wv = GetGoodFFTSizes()
+	numEntries = DimSize(wv, ROWS)
+	CHECK_GT_VAR(numEntries, 0)
+
+	Make/FREE/N=(numEntries) junkWave
+
+	junkWave[] = CheckPrimeFactors(wv[p])
+End