diff --git a/src/classes/PFourier.cpp b/src/classes/PFourier.cpp index 924a6f956..0beb67de7 100644 --- a/src/classes/PFourier.cpp +++ b/src/classes/PFourier.cpp @@ -51,7 +51,18 @@ // Constructor //-------------------------------------------------------------------------- /** + *

Default constructor for phase correction optimizer. * + *

Initializes a phase correction object with default Fourier data. + * The actual Fourier data must be set externally before minimization. + * + * @param minBin Minimum frequency bin index for optimization region. + * Use -1 to optimize over all bins. Default: -1. + * @param maxBin Maximum frequency bin index for optimization region. + * Use -1 to optimize over all bins. Default: -1. + * + *

Note: Restricting the optimization range to signal-containing + * regions improves results by excluding noisy baseline regions. */ PFTPhaseCorrection::PFTPhaseCorrection(const Int_t minBin, const Int_t maxBin) : fMinBin(minBin), fMaxBin(maxBin) @@ -63,7 +74,28 @@ PFTPhaseCorrection::PFTPhaseCorrection(const Int_t minBin, const Int_t maxBin) : // Constructor //-------------------------------------------------------------------------- /** + *

Constructor with explicit Fourier transform data. * + *

Creates a phase correction optimizer with pre-computed Fourier data. + * The optimizer will find phase parameters φ(ω) = c₀ + c₁·ω that maximize + * the real component while minimizing the imaginary component. + * + * @param reFT Real part of the Fourier transform (vector of amplitudes). + * @param imFT Imaginary part of the Fourier transform (vector of amplitudes). + * @param minBin Minimum frequency bin index for optimization region. + * Use -1 to optimize over all bins. Default: -1. + * @param maxBin Maximum frequency bin index for optimization region. + * Use -1 to optimize over all bins. Default: -1. + * + *

Usage example: + * @code + * std::vector real, imag; + * // ... fill real and imag with FFT results ... + * PFTPhaseCorrection phCorr(real, imag, 10, 500); + * phCorr.Minimize(); + * Double_t c0 = phCorr.GetPhaseCorrectionParam(0); + * Double_t c1 = phCorr.GetPhaseCorrectionParam(1); + * @endcode */ PFTPhaseCorrection::PFTPhaseCorrection(std::vector &reFT, std::vector &imFT, const Int_t minBin, const Int_t maxBin) : fReal(reFT), fImag(imFT), fMinBin(minBin), fMaxBin(maxBin) @@ -86,7 +118,28 @@ PFTPhaseCorrection::PFTPhaseCorrection(std::vector &reFT, std::vector< // Minimize (public) //-------------------------------------------------------------------------- /** + *

Performs phase correction minimization using Minuit2. * + *

This method finds optimal phase parameters (c₀, c₁) by minimizing + * a combined entropy-penalty functional: + * - Entropy term: Promotes smooth, concentrated real spectra + * - Penalty term: Penalizes negative values in real spectrum + * + *

The phase correction is: φ(ω) = c₀ + c₁·ω where: + * - c₀ corrects constant phase offset (time-zero uncertainty) + * - c₁ corrects linear phase dispersion (detector timing effects) + * + *

Algorithm: + * 1. Initialize Minuit2 parameters with starting values (c₀=0, c₁=0) + * 2. Minimize entropy + γ×penalty using MnMinimize + * 3. Store optimal parameters in fPh_c0, fPh_c1 + * + *

After calling this method, use GetPhaseCorrectionParam() to retrieve + * the optimal phase parameters and IsValid() to check success. + * + * @see GetPhaseCorrectionParam() + * @see SetGamma() + * @see IsValid() */ void PFTPhaseCorrection::Minimize() { @@ -118,7 +171,26 @@ void PFTPhaseCorrection::Minimize() // GetPhaseCorrectionParam (public) //-------------------------------------------------------------------------- /** + *

Retrieves optimized phase correction parameter by index. * + *

Returns the phase parameter determined by Minimize(). The phase + * correction formula is: φ(ω) = c₀ + c₁·(ω/ω_max) + * + * @param idx Parameter index: + * - 0: Returns c₀ (constant phase offset in radians) + * - 1: Returns c₁ (linear phase dispersion coefficient) + * + * @return Phase correction parameter value, or 0.0 if idx is invalid. + * + *

Example: + * @code + * phCorr.Minimize(); + * if (phCorr.IsValid()) { + * Double_t c0 = phCorr.GetPhaseCorrectionParam(0); + * Double_t c1 = phCorr.GetPhaseCorrectionParam(1); + * std::cout << "Phase: " << c0 << " + " << c1 << "*w" << std::endl; + * } + * @endcode */ Double_t PFTPhaseCorrection::GetPhaseCorrectionParam(UInt_t idx) { @@ -138,7 +210,15 @@ Double_t PFTPhaseCorrection::GetPhaseCorrectionParam(UInt_t idx) // GetMinimum (public) //-------------------------------------------------------------------------- /** + *

Returns the minimum value of the optimization functional. * + *

This value represents the combined entropy-penalty functional at + * the optimal phase parameters. Lower values indicate better phase + * correction. A value of -1.0 indicates failed minimization. + * + * @return Minimum functional value, or -1.0 if minimization failed. + * + *

Note: Always check IsValid() before using this value. */ Double_t PFTPhaseCorrection::GetMinimum() { @@ -154,7 +234,16 @@ Double_t PFTPhaseCorrection::GetMinimum() // Init (private) //-------------------------------------------------------------------------- /** + *

Initializes phase correction object with default values. * + *

Sets initial state: + * - fValid = true (object ready for use) + * - fPh_c0 = 0.0 (no constant phase offset) + * - fPh_c1 = 0.0 (no linear phase dispersion) + * - fGamma = 1.0 (equal weighting of entropy and penalty) + * - fMin = -1.0 (no minimization performed yet) + * + *

Called by constructors to establish consistent initial state. */ void PFTPhaseCorrection::Init() { @@ -169,7 +258,20 @@ void PFTPhaseCorrection::Init() // CalcPhasedFT (private) //-------------------------------------------------------------------------- /** + *

Calculates phase-corrected Fourier transform. * + *

Applies the phase correction φ(ω) = c₀ + c₁·(i/N) to the complex + * Fourier spectrum using rotation in the complex plane: + * - F'_re(ω) = F_re(ω)·cos(φ) - F_im(ω)·sin(φ) + * - F'_im(ω) = F_re(ω)·sin(φ) + F_im(ω)·cos(φ) + * + *

The corrected spectra are stored in fRealPh and fImagPh. + * + *

Physics: This rotation compensates for: + * - Time-zero offset (constant phase c₀) + * - Timing dispersion (linear phase c₁) + * + * @see CalcRealPhFTDerivative() */ void PFTPhaseCorrection::CalcPhasedFT() const { @@ -188,7 +290,19 @@ void PFTPhaseCorrection::CalcPhasedFT() const // CalcRealPhFTDerivative (private) //-------------------------------------------------------------------------- /** + *

Calculates first derivative of phase-corrected real Fourier spectrum. * + *

Computes the finite difference derivative: dF/dω ≈ F(i+1) - F(i) + * for all interior points. Boundary points are set to 1.0. + * + *

The derivative is used in the entropy calculation, where smooth + * spectra (small derivatives) have lower entropy and are preferred by + * the optimization algorithm. + * + *

Purpose: Entropy minimization favors phase corrections that + * produce smooth, well-resolved real spectra with minimal oscillations. + * + * @see Entropy() */ void PFTPhaseCorrection::CalcRealPhFTDerivative() const { @@ -205,7 +319,22 @@ void PFTPhaseCorrection::CalcRealPhFTDerivative() const // Penalty (private) //-------------------------------------------------------------------------- /** + *

Calculates penalty term for negative real spectrum values. * + *

Computes: P = γ·Σ[F_re(ω)² for all F_re(ω) < 0] + * + *

The penalty term penalizes negative values in the real Fourier + * spectrum, since physically meaningful absorption-mode spectra should + * be predominantly positive. The γ parameter controls the relative + * importance of this constraint. + * + * @return Weighted penalty value (γ×sum of squared negative amplitudes) + * + *

Note: Only bins in range [fMinBin, fMaxBin) are considered, + * allowing focus on signal-containing regions. + * + * @see SetGamma() + * @see Entropy() */ Double_t PFTPhaseCorrection::Penalty() const { @@ -223,7 +352,26 @@ Double_t PFTPhaseCorrection::Penalty() const // Entropy (private) //-------------------------------------------------------------------------- /** + *

Calculates Shannon entropy of the real spectrum derivative. * + *

Computes: S = -Σ[p_i·ln(p_i)] where p_i = |dF/dω|_i / Σ|dF/dω| + * + *

This entropy measure quantifies the "smoothness" of the real + * Fourier spectrum: + * - Low entropy: Concentrated, smooth spectrum (desired) + * - High entropy: Dispersed, noisy spectrum (undesired) + * + *

Physical interpretation: Correct phase alignment produces + * sharp, well-defined peaks with smooth baselines, resulting in + * concentrated derivatives and low entropy. + * + * @return Shannon entropy of normalized derivative distribution + * + *

Note: Small values (< 10⁻¹⁵) are skipped to avoid + * numerical issues with log(0). + * + * @see CalcRealPhFTDerivative() + * @see Penalty() */ Double_t PFTPhaseCorrection::Entropy() const { @@ -250,7 +398,29 @@ Double_t PFTPhaseCorrection::Entropy() const // operator() (private) //-------------------------------------------------------------------------- /** + *

Objective function for Minuit2 minimization (FCNBase interface). * + *

Evaluates the combined entropy-penalty functional: + * f(c₀, c₁) = S + γ·P + * where: + * - S = entropy of phase-corrected real spectrum derivative + * - P = penalty for negative values in real spectrum + * - γ = balancing parameter (fGamma) + * + * @param par Parameter vector: [0]=c₀ (constant phase), [1]=c₁ (linear phase) + * @return Combined functional value to be minimized + * + *

Algorithm flow: + * 1. Apply phase correction with parameters par + * 2. Calculate phase-corrected Fourier transform + * 3. Compute spectrum derivative + * 4. Return entropy + penalty + * + *

Called repeatedly by Minuit2 during optimization. + * + * @see Entropy() + * @see Penalty() + * @see Minimize() */ double PFTPhaseCorrection::operator()(const std::vector &par) const { @@ -271,15 +441,43 @@ double PFTPhaseCorrection::operator()(const std::vector &par) const // Constructor //-------------------------------------------------------------------------- /** - *

Constructor. + *

Constructor for Fourier transform engine. * - * \param data data histogram - * \param unitTag tag telling in which units the Fourier transform shall be represented. Possible tags are: - * FOURIER_UNIT_GAUSS, FOURIER_UNIT_TESLA, FOURIER_UNIT_FREQ, FOURIER_UNIT_CYCLES - * \param startTime start time of the data time window - * \param endTime end time of the data time window - * \param dcCorrected if true, removed DC offset from signal before Fourier transformation, otherwise not - * \param zeroPaddingPower if set to values > 0, there will be zero padding up to 2^zeroPaddingPower + *

Creates a PFourier object configured to transform μSR time-domain + * data to frequency/field domain. The constructor validates inputs, + * allocates FFTW resources, and prepares for FFT computation. + * + * @param data TH1F histogram containing time-domain μSR signal. + * X-axis: time in microseconds, Y-axis: asymmetry or counts. + * @param unitTag Output units for frequency axis: + * - FOURIER_UNIT_GAUSS (1): Magnetic field in Gauss + * - FOURIER_UNIT_TESLA (2): Magnetic field in Tesla + * - FOURIER_UNIT_FREQ (3): Frequency in MHz + * - FOURIER_UNIT_CYCLES (4): Angular frequency in Mc/s (Mrad/s) + * @param startTime Start of transform window in μs (0.0 = from t₀). Default: 0.0 + * @param endTime End of transform window in μs (0.0 = to end). Default: 0.0 + * @param dcCorrected If true, remove DC offset before FFT to eliminate + * zero-frequency artifact. Default: false + * @param zeroPaddingPower Zero-pad to 2^N points where N=zeroPaddingPower. + * Value 0 disables padding. Padding improves frequency + * interpolation but doesn't add information. Default: 0 + * + *

Unit conversions (using γ_μ/2π = 0.01355 MHz/G): + * - Gauss: B(G) = ω(MHz) / 0.01355 + * - Tesla: B(T) = ω(MHz) / 135.54 + * + *

Example usage: + * @code + * TH1F *timeHisto = ...; // μSR time spectrum + * PFourier ft(timeHisto, FOURIER_UNIT_GAUSS, 0.1, 10.0, true, 12); + * ft.Transform(F_APODIZATION_WEAK); + * TH1F *powerSpec = ft.GetPowerFourier(); + * @endcode + * + *

After construction, check IsValid() before calling Transform(). + * + * @see Transform() + * @see IsValid() */ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTime, Bool_t dcCorrected, UInt_t zeroPaddingPower, Bool_t useFFTW) : fData(data), fUnitTag(unitTag), fStartTime(startTime), fEndTime(endTime), @@ -429,7 +627,16 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi // Destructor //-------------------------------------------------------------------------- /** - *

Destructor + *

Destructor - cleans up FFTW resources. + * + *

Releases all FFTW-allocated memory and destroys the FFT plan. + * Proper cleanup is essential to avoid memory leaks, as FFTW uses + * special aligned memory allocation. + * + *

Resources freed: + * - fFFTwPlan: FFTW execution plan + * - fIn: Input array (time-domain data) + * - fOut: Output array (frequency-domain data) */ PFourier::~PFourier() { @@ -458,10 +665,37 @@ PFourier::~PFourier() // Transform (public) //-------------------------------------------------------------------------- /** - *

Carries out the Fourier transform. It is assumed that fStartTime is the time zero - * for the Fourier frame. Hence if fStartTime != 0.0 the phase shift will be corrected. + *

Executes the Fourier transform with optional apodization. * - * \param apodizationTag 0=no apod., 1=weak apod., 2=medium apod., 3=strong apod. + *

This is the main workhorse method that performs the complete FFT + * pipeline: + * 1. Prepare input data (DC correction, zero padding) + * 2. Apply apodization window if requested + * 3. Execute FFTW plan (compute FFT) + * 4. Correct phase for non-zero start time + * + *

Phase correction: If fStartTime ≠ 0, applies phase shift + * e^(i·ω·t₀) to account for the time-zero offset. This ensures proper + * phase relationships in the frequency spectrum. + * + *

Apodization: Window functions reduce spectral leakage + * (Gibbs phenomenon) at the cost of slight peak broadening: + * - NONE: Best amplitude accuracy, worst leakage + * - WEAK: Minimal smoothing + * - MEDIUM: Balanced (recommended for most cases) + * - STRONG: Maximum leakage suppression + * + * @param apodizationTag Window function strength: + * - 0 or F_APODIZATION_NONE: No windowing + * - 1 or F_APODIZATION_WEAK: Gentle taper + * - 2 or F_APODIZATION_MEDIUM: Moderate taper + * - 3 or F_APODIZATION_STRONG: Heavy taper + * + *

After calling Transform(), retrieve results using: + * GetRealFourier(), GetImaginaryFourier(), GetPowerFourier(), GetPhaseFourier() + * + * @see PrepareFFTwInputData() + * @see GetPowerFourier() */ void PFourier::Transform(UInt_t apodizationTag) { @@ -542,7 +776,25 @@ void PFourier::SetUseFFTW(const Bool_t flag) // GetMaxFreq (public) //-------------------------------------------------------------------------- /** - *

returns the maximal frequency in units choosen, i.e. Gauss, Tesla, MHz, Mc/s + *

Returns the maximum frequency (Nyquist frequency) of the spectrum. + * + *

The Nyquist frequency is the highest frequency that can be + * unambiguously represented given the time resolution: + * f_Nyquist = 1 / (2·Δt) + * + *

Higher frequencies are aliased back into the [0, f_Nyquist] range, + * so the useful spectrum extends from 0 to f_Nyquist. + * + * @return Maximum frequency in units specified by fUnitTag: + * - Gauss (if FOURIER_UNIT_GAUSS) + * - Tesla (if FOURIER_UNIT_TESLA) + * - MHz (if FOURIER_UNIT_FREQ) + * - Mc/s (if FOURIER_UNIT_CYCLES) + * + *

Example: For time resolution Δt = 0.01 μs: + * f_Nyquist = 1/(2×0.01) = 50 MHz ≈ 3690 Gauss + * + * @see GetResolution() */ Double_t PFourier::GetMaxFreq() { @@ -555,9 +807,36 @@ Double_t PFourier::GetMaxFreq() // GetRealFourier (public) //-------------------------------------------------------------------------- /** - *

returns the real part Fourier as a histogram. + *

Returns real part of Fourier transform as ROOT histogram. * - * \param scale normalisation factor + *

The real component represents the in-phase (cosine) projection + * of the signal. Without phase correction, both positive and negative + * values typically appear due to arbitrary phase offsets. + * + *

Interpretation: + * - Phase-corrected: Absorption-mode spectrum (predominantly positive) + * - Uncorrected: Mixed-phase spectrum (positive + negative regions) + * + *

Histogram covers [0, f_Nyquist] with N/2 bins where N is the + * FFT size (including zero padding if applied). + * + * @param scale Multiplication factor applied to all amplitudes. + * Use for normalization or unit conversion. Default: 1.0 + * + * @return Pointer to newly allocated TH1F histogram. Caller owns + * the histogram and must delete it. Returns nullptr on error. + * + *

Example: + * @code + * ft.Transform(F_APODIZATION_WEAK); + * TH1F *realSpec = ft.GetRealFourier(1.0); + * realSpec->Draw(); + * // ... use histogram ... + * delete realSpec; + * @endcode + * + * @see GetImaginaryFourier() + * @see GetPhaseOptRealFourier() */ TH1F* PFourier::GetRealFourier(const Double_t scale) { @@ -602,16 +881,54 @@ TH1F* PFourier::GetRealFourier(const Double_t scale) // GetPhaseOptRealFourier (public, static) //-------------------------------------------------------------------------- /** - *

returns the phase corrected real Fourier transform. + *

Returns phase-optimized real Fourier spectrum (static method). * - * \return the TH1F histo of the phase 'optimzed' real Fourier transform. + *

This static method performs automatic phase correction on a complex + * Fourier spectrum to maximize the real component (absorption mode) while + * minimizing the imaginary component. The algorithm finds optimal phase + * parameters φ(ω) = c₀ + c₁·ω that produce the cleanest real spectrum. * - * \param re real part Fourier histogram - * \param im imaginary part Fourier histogram - * \param phase return value of the optimal phase dispersion phase[0]+phase[1]*i/N - * \param scale normalisation factor - * \param min minimal freq / field from which to optimise. Given in the choosen unit. - * \param max maximal freq / field up to which to optimise. Given in the choosen unit. + *

Why phase correction? Raw FFT produces mixed-phase spectra + * with arbitrary phase offsets due to: + * - Uncertain time-zero (t₀) determination + * - Detector timing offsets + * - Signal propagation delays + * + *

Algorithm: + * 1. Extract Re/Im data in specified frequency range [min, max] + * 2. Use PFTPhaseCorrection to minimize entropy+penalty functional + * 3. Apply optimal phase rotation to entire spectrum + * 4. Return phase-corrected real part + * + *

Optimization range [min, max]: Restrict to signal-containing + * regions for best results. Including noisy baseline degrades optimization. + * + * @param re Real part of Fourier transform (TH1F histogram) + * @param im Imaginary part of Fourier transform (TH1F histogram) + * @param phase Output parameter: phase[0]=c₀, phase[1]=c₁ (vector resized to 2) + * @param scale Amplitude scaling factor. Default: 1.0 + * @param min Minimum frequency/field for optimization window. + * Use -1.0 to include all frequencies. Default: -1.0 + * @param max Maximum frequency/field for optimization window. + * Use -1.0 to include all frequencies. Default: -1.0 + * + * @return Pointer to newly allocated phase-corrected TH1F histogram. + * Caller owns and must delete. Returns nullptr on error. + * + *

Example: + * @code + * TH1F *real = ft.GetRealFourier(); + * TH1F *imag = ft.GetImaginaryFourier(); + * std::vector phaseParams; + * TH1F *phCorrected = PFourier::GetPhaseOptRealFourier(real, imag, phaseParams, 1.0, 5.0, 50.0); + * std::cout << "Phase: " << phaseParams[0] << " + " << phaseParams[1] << "*w" << std::endl; + * phCorrected->Draw(); + * delete real; delete imag; delete phCorrected; + * @endcode + * + * @see PFTPhaseCorrection + * @see GetRealFourier() + * @see GetImaginaryFourier() */ TH1F* PFourier::GetPhaseOptRealFourier(const TH1F *re, const TH1F *im, std::vector &phase, const Double_t scale, const Double_t min, const Double_t max) @@ -703,9 +1020,29 @@ TH1F* PFourier::GetPhaseOptRealFourier(const TH1F *re, const TH1F *im, std::vect // GetImaginaryFourier (public) //-------------------------------------------------------------------------- /** - *

returns the imaginary part Fourier as a histogram. + *

Returns imaginary part of Fourier transform as ROOT histogram. * - * \param scale normalisation factor + *

The imaginary component represents the out-of-phase (sine) + * projection of the signal. Ideally zero for perfect phase correction, + * but contains signal information when phase is uncorrected. + * + *

Uses: + * - Phase correction optimization (minimize imaginary component) + * - Dispersion-mode spectroscopy + * - Computing power spectrum: |F|² = Re² + Im² + * - Phase spectrum calculation: φ = atan(Im/Re) + * + *

Histogram covers [0, f_Nyquist] with N/2 bins where N is the + * FFT size (including zero padding if applied). + * + * @param scale Multiplication factor applied to all amplitudes. + * Use for normalization or unit conversion. Default: 1.0 + * + * @return Pointer to newly allocated TH1F histogram. Caller owns + * the histogram and must delete it. Returns nullptr on error. + * + * @see GetRealFourier() + * @see GetPhaseOptRealFourier() */ TH1F* PFourier::GetImaginaryFourier(const Double_t scale) { @@ -750,9 +1087,37 @@ TH1F* PFourier::GetImaginaryFourier(const Double_t scale) // GetPowerFourier (public) //-------------------------------------------------------------------------- /** - *

returns the Fourier power spectrum as a histogram. + *

Returns power spectrum |F(ω)| as ROOT histogram. * - * \param scale normalisation factor + *

The power spectrum shows signal magnitude at each frequency, + * computed as: |F(ω)| = √(Re² + Im²) + * + *

Advantages: + * - Always positive (easier interpretation) + * - Phase-independent (no phase correction needed) + * - Directly shows signal strength vs. frequency + * - Ideal for identifying all frequency components + * + *

Applications: + * - Quick survey of frequency content + * - Field distribution analysis in superconductors + * - Multiple muon site identification + * - TF-μSR field measurements + * + *

Histogram covers [0, f_Nyquist] with N/2 bins where N is the + * FFT size (including zero padding if applied). + * + * @param scale Multiplication factor applied to all amplitudes. + * Use for normalization or unit conversion. Default: 1.0 + * + * @return Pointer to newly allocated TH1F histogram. Caller owns + * the histogram and must delete it. Returns nullptr on error. + * + *

Note: This returns the amplitude |F|, not the power |F|². + * For power values, square the histogram contents. + * + * @see GetRealFourier() + * @see GetPhaseFourier() */ TH1F* PFourier::GetPowerFourier(const Double_t scale) { @@ -797,9 +1162,39 @@ TH1F* PFourier::GetPowerFourier(const Double_t scale) // GetPhaseFourier (public) //-------------------------------------------------------------------------- /** - *

returns the Fourier phase spectrum as a histogram. + *

Returns phase spectrum φ(ω) = arg(F) as ROOT histogram. * - * \param scale normalisation factor + *

The phase spectrum shows the phase angle of the complex Fourier + * transform at each frequency: φ(ω) = atan2(Im, Re) + * + *

Range: Phase values are in [-π, +π] radians, with proper + * quadrant handling using the two-argument arctangent. + * + *

Interpretation: + * - Constant phase: All frequencies in phase (simple exponential decay) + * - Linear phase: Time-zero offset (shift theorem) + * - Quadratic phase: Dispersion or chirp + * - Random phase: Noise + * + *

Applications: + * - Diagnosing time-zero problems + * - Verifying phase correction quality + * - Identifying signal vs. noise regions + * + *

Histogram covers [0, f_Nyquist] with N/2 bins where N is the + * FFT size (including zero padding if applied). + * + * @param scale Multiplication factor applied to phase values. + * Typically 1.0 (radians) or 180/π (degrees). Default: 1.0 + * + * @return Pointer to newly allocated TH1F histogram. Caller owns + * the histogram and must delete it. Returns nullptr on error. + * + *

Note: Phase is undefined where amplitude is zero. In practice, + * phase values in noise regions are meaningless. + * + * @see GetRealFourier() + * @see GetImaginaryFourier() */ TH1F* PFourier::GetPhaseFourier(const Double_t scale) { @@ -865,11 +1260,32 @@ TH1F* PFourier::GetPhaseFourier(const Double_t scale) // PrepareFFTwInputData (private) //-------------------------------------------------------------------------- /** - *

Feeds the Fourier data and apply the apodization. + *

Prepares time-domain data for FFT with optional DC correction and apodization. * - * \param apodizationTag apodization tag. Possible are currently: F_APODIZATION_NONE = no apodization, - * F_APODIZATION_WEAK = weak apodization, F_APODIZATION_MEDIUM = intermediate apodization, - * F_APODIZATION_STRONG = strong apodization + *

This method performs essential preprocessing steps: + * 1. Locates time-zero bin in the histogram + * 2. Extracts data from [fStartTime, fEndTime] window + * 3. Removes DC offset (if fDCCorrected=true) for baseline correction + * 4. Zero-pads to fNoOfBins (power of 2 for optimal FFT performance) + * 5. Applies apodization window to reduce spectral leakage + * + *

DC Correction: Subtracts mean value from time signal to + * remove constant offset, which would otherwise create a large artifact + * at zero frequency. + * + *

Zero Padding: Increases frequency resolution by interpolation + * without adding information content. Often used for smoother spectra. + * + * @param apodizationTag Window function strength: + * - F_APODIZATION_NONE (1): Rectangular window (no smoothing) + * - F_APODIZATION_WEAK (2): Gentle roll-off at edges + * - F_APODIZATION_MEDIUM (3): Moderate smoothing + * - F_APODIZATION_STRONG (4): Heavy smoothing for maximum leakage suppression + * + *

Result is stored in fIn[] array ready for FFTW execution. + * + * @see ApodizeData() + * @see Transform() */ void PFourier::PrepareFFTwInputData(UInt_t apodizationTag) { @@ -928,11 +1344,33 @@ void PFourier::PrepareFFTwInputData(UInt_t apodizationTag) // ApodizeData (private) //-------------------------------------------------------------------------- /** - *

Carries out the appodization of the data. + *

Applies apodization (windowing) to time-domain data before FFT. * - * \param apodizationTag apodization tag. Possible are currently: F_APODIZATION_NONE = no apodization, - * F_APODIZATION_WEAK = weak apodization, F_APODIZATION_MEDIUM = intermediate apodization, - * F_APODIZATION_STRONG = strong apodization + *

Purpose: Apodization multiplies the time signal by a smooth + * window function that tapers to zero at the edges, reducing the Gibbs + * phenomenon (spectral leakage) caused by finite observation windows. + * + *

Trade-off: + * - Advantage: Sharper, cleaner frequency peaks with less sidelobes + * - Disadvantage: Slightly broader main peaks, reduced amplitude accuracy + * + *

Window functions: Polynomial windows of the form: + * w(t) = Σ c_j·(t/T)^(2j) where t ∈ [0,T] + * + *

Three predefined window strengths with optimized coefficients: + * - Weak: Minimal smoothing, preserves amplitude + * - Medium: Balanced smoothing for general use + * - Strong: Maximum leakage suppression for crowded spectra + * + * @param apodizationTag Window function type: + * - F_APODIZATION_NONE: No windowing (rectangular) + * - F_APODIZATION_WEAK: Gentle taper + * - F_APODIZATION_MEDIUM: Moderate taper + * - F_APODIZATION_STRONG: Heavy taper + * + *

Window is applied in-place to fIn[] array. + * + * @see PrepareFFTwInputData() */ void PFourier::ApodizeData(Int_t apodizationTag) {