function spOpt = spOptions(avsp)
% This script sets options and parameters for signal processing performed
% by "mainTOFsPROc.m"
%
% OUTPUT:
% spOpt - structure with option and parameter settings in the fields
%
% INPUT:
% avsp - average spectrum for the studied data set
%
% USAGE:
% spopt=spOptions(avsp);
% Dependency: none
% NOTE: Current default settings are optimized for the broad-range linear
% Ultraflex II MALDI-TOF spectra [Bunai et al: JPR 2007],
% e.g. included in "qc11dat.mat": "spmqc11".
%
% For new TOF data, the parameters should be optimized using average
% spectrum for the data set and following instructions for each field
% below using corresponding auxiliary functions of the toolbox or built-ins.
% For input/output and usage description for auxiliary functions, type:
% help function_name;
%
%% GLOBALLY USED PARAMETERS
% mass spectrometer ADC precision
spOpt.precn =8; % 8 - bit ADC for Ultraflex and PBS instruments
%% BASELINE SUBTRACTION PARAMETERS
% choice of baseline model: 1 - charge accumulation; 0 (default) - analytical
% function (Gaussian, exponential or linear)
spOpt.modBS = 0;
% charge accumulation model is appropriate if long decay tails are observed
% for peaks (in the presence of coupling of detector amplifier)
% see [Malyarenko et al: Clin Chem 2005]
% choice of analytical function
spOpt.chBM = 2; % 1 - linear, 2 - exponential, 3 - Gaussian
% to determine appropriate choice: plot(log(abs(avsp(first_half)-offset))
% offset = mean(avsp(spOpt.osS)); % constant offset
% exponential baseline would produce linear log-trend, and Gaussian -
% parabolic in early data range
% parameters for charge accumulation model (spOpt.modBS=1)
% see [Malyarenko et al: Clin Chem 2005]
spOpt.amp=0.0007; % amplitude factor
spOpt.decay = 0.9989; % decay rate
% try applying to the average spectrum using "bscharge.m" and observe
% correction of the slowly decaying right peak side (if appropriate)
% TOF range to estimate constant ADC offset (late in the record: far from
% a matrix hump and without peaks)
spOpt.osS = 22300:22800; % is used to estimate the offset = mean(osS)
% range of spectrum points without peaks on the decaying part (after the hump,
% before the leveling off) for baseline-trend estimate by "bsModel.m".
% estimate by visual inspection: plot(avsp);
spOpt.bsS=4; % use "auto-trend" %[5370:5420,6600:6700];
% NOTE: *automated* trend estimate can be attempted, when set to a single number (length(bsS)=1)
% to check baseline subtraction performance:
% plot([avsp, baseline]); avsp_bs = avsp - baseline;
% NOTE: when pedestal removal is used downstream, the baseline model errors
% are not significant for the overall performance
%% INTEGRATIVE DOWN-SAMPLING PARAMETERS
% choose high-frequency noise model to determine if it is necessary to scale
% the down-sampled noise by sqrt(window)
spOpt.noiseMod = 0; % 1 - Gaussian (white, random) noise; 0 - derivative noise (also default)
% derivative noise (e.g., PBS: produced by differential amplifier) exhibits
% higher FFT spectral density for lower frequencies; spectral density of
% white noise is uniform. Use Gaussian, when not sure.
% choose minimal peak width desired after down-sampling (in time clicks)
spOpt.minHW = 4.5; % minimal acceptable peak width: 3-5 is recommended
% constant peak width observed over the mass focusing range (in time clicks)
spOpt.cHW = 5;
% find left-half width of a peak in early mass range, e.g., by visual
% inspection: plot(avsp-baseline);
% or using "peakwidth_pf.m" output in the constant peak-width range
% parameters for quadratic fit of peak half-width dependence on TOF
spOpt.pfit =[3.9338e-7,-0.0066998,33.278]; % quadratic polynomial coefficients: [p1,p2,p3]
% apply "peakwidth_pf.m" for (avsp-baseline) to find observed peak
% half-width dependence on time; remove the outliers (usually, of higher
% width) due to multiple overlapping peaks or noise; and perform a
% quadratic fit to find the coefficients:
% polyfit(peak_positions, peak_width, 2); t = (1:length(avsp))';
% an example of the "peakwidth_pf.m" output: [peak_position, peak_width]
% pair for the deafult paramepetrs is provide with this toolbox in "hwqc11" array
% (saved in working/dataexmpl/qc11dat_spm_rmz_hw.mat)
% width_fit = p1*t.^2+p2*t+p3;
% to check performance visually:
% plot(t, [avsp*max(peak_width)/max(avsp), width_fit], peak_positions, peak_width, 'rs');
% NOTE: to check the quality of integrative down-sampling, after running
%"resample_rev2.m", apply "peakwidth_pf.m" to the down sampled spectrum:
% the result should produce a constant peak-width (pHW) over the full data range.
% Optimal performance may require a couple iterations through the "pfit" parameters
% and performance check after down-sampling.
%% FILTER PARAMETERS AND OPTIONS
% choose the filter: see [Malyarenko et al: RCMS 2006] for guidelines
spOpt.chFilter = 1; % 0 - none; 1 - OLF (also default); 2 - non-linear; 3 - matched; 4 - MAV
% NOTE: assumes that peak-width is constant (down-sampling is required before
% application in the broad m/z range)
% peak wavelet model choice for target filter (OLF or non-linear):
% left-Gaussian-right-Lorentzian (GL), asymmetric Gaussian (GG), or asymmetric Lorentzian (LL)
spOpt.peakMdl = 2; % 1 - GL wavelet (also default); 2 - GG wavelet; 3- LL wavelet
% to find best analytical approximation for the peak shape, choose a single
% peak in the spectrum: mpk = avsp_ds(peak_range); mpk=mpk/max(mpk); % to normalize
% and fit it to the analytical shapes: GL, GG or LL, e.g., using built-in "normpdf.m" to
% model a Gaussian, and "lorf.m" for Lorentzian, e.g.: wvlt =
% [gaussian(1:max_position)/max(gaussian), lorentzian(max_position+1:end)/max_lorentzian];
% choose the best fit by minimizing the difference between experiment and
% model: plot(mpk-wvlt);
% NOTE: inappropriate peak shape may cause mis-shaping artifacts by target filtering
% by separately fitting the left and right half of the model peak,
% determine if the right-side has different width:
% lasym = right_sigma -left_sigma
spOpt.lasym = 0;% right-side asymmetry for a peak model
% can be negative, if left-half is wider
%% PEDESTAL REMOVAL PARAMETERS
% NOTE: assumes constant peak-width, should be used on down-sampled data (after IDS)
spOpt.nfPed = 10; % number fraction of peak half-width for the MAV window of local minima trend
% higher fractions will produce smoother trends (less of peak pedestals
% will be removed). Inspect visually for a few choices: 5, 10, 20 and
% with above estimated half-width of a peak (pHW) in down-sampled/filtered signal:
% spd = bsMinMav(avsp_ids, pHW, 10); plot ([avsp, spd]); avsp_pr = avsp-spd;
%% PEAK DETECTION PARAMETERS
% choose the trustworthy range, above matrix or mass deflector artifacts, or residual ringing
% from detector overload by visual inspection of the average spectrum after pedestal removal
% plot(avsp_pr);
% can be estimated as (original deflector/matrix hump time)/(constant down sampling rate)
spOpt.startPeak = 5900; % onset of peak detection in down-sampled time ticks (TOF)
% visually check if high-frequency noise amplitude changes with TOF: plot(avsp);
% can be evaluated from: std(avsp_bs(range)); % using different opt.bsS
% ranges without peaks, and modeled analytically: see "noise amplitude" in "mainTOFsPROc. m"
spOpt.nsdepBS = 0; % set to 1 if noise is sqrt-dependent (Poisson) on baseline amplitude
% e.g., use 1 for Bruker Ultraflex TOF, and 0 for Ciphergen PBS TOF
% used for estimate of noise amplitude after different filtering options
% if not sure, use 0 for stationary noise (constant noise amplitude) -- may
% produce higher false detection in the region of higher noise amplitude,
% and require manual correction of the peak list
% Signal-to-noise (amplitude) threshold for peak detection in respect to global
% baseline and local minima, used by "trivialPeakFinderSNR.m"
% noise amplitude parameter is estimated automatically depending on
% spOpt.nsdepBS and filter choice: see "noise amplitude" in "mainTOFsPROc. m"
% to visually check the quality of peak detection (false and true peaks)
% plot(1:length(avsp_pr), avsp_pr, pks, avsp_pr(pks), 'ko');
spOpt.SNR = 4;% SNR threshold
