gavindiaz/QuanTAlib
1
0
Fork 0
mirror of https://github.com/mihakralj/QuanTAlib.git synced 2026-08-22 20:48:04 +00:00
Code Issues Packages Projects Releases Wiki Activity

fix(docs): correct .md documentation across errors, dynamics, filters, forecasts, momentum, numerics, oscillators, reversals, statistics, trends, volatility, volume

Browse Source 
Deep review of all indicator categories verified .md headers against .cs WarmupPeriod, parameters, inputs, and outputs. Fixes include warmup corrections, parameter documentation, output type accuracy, and Pine Script alignment.
This commit is contained in:
Miha Kralj
2026-03-10 18:38:23 -07:00
parent 8906c62dcf
commit 35a6702b06
178 changed files with 2579 additions and 998 deletions
Download Patch File Download Diff File Expand all files Collapse all files
lib/numerics/fft/Fft.cs
+232 -109
View File
@@ -1,8 +1,11 @@
// FFT: Fast Fourier Transform — Dominant Cycle Detector
// Estimates the dominant cycle period in bars using a DFT on a windowed price buffer.
// Algorithm: Ehlers, J.F. "Cycle Analytics for Traders." Wiley, 2013.
// Hanning-windowed DFT across bins [minBin..maxBin], with parabolic interpolation
// for sub-bin period estimation. Output: dominant cycle period in bars (clamped).
// Estimates the dominant cycle period in bars using a radix-2 Cooley-Tukey FFT
// on a Hanning-windowed price buffer, with parabolic interpolation for sub-bin
// period estimation. Output: dominant cycle period in bars (clamped).
//
// Algorithm: Cooley, J.W. & Tukey, J.W. (1965). "An Algorithm for the Machine
// Calculation of Complex Fourier Series." Mathematics of Computation, 19(90).
// Ehlers, J.F. "Cycle Analytics for Traders." Wiley, 2013 (application context).
using System.Buffers;
using System.Runtime.CompilerServices;
@@ -12,15 +15,16 @@ namespace QuanTAlib;
/// <summary>
/// FFT: Fast Fourier Transform Dominant Cycle Detector
/// Computes the dominant cycle period using a Hanning-windowed DFT
/// over a rolling price buffer, with parabolic interpolation refinement.
/// Computes the dominant cycle period using a Hanning-windowed radix-2
/// Cooley-Tukey FFT over a rolling price buffer, with parabolic interpolation.
/// </summary>
/// <remarks>
/// Key properties:
/// - Output: dominant cycle period in bars, clamped to [minPeriod, maxPeriod]
/// - windowSize must be 32, 64, or 128
/// - windowSize must be 32, 64, or 128 (power of 2 for radix-2)
/// - WarmupPeriod = windowSize bars
/// - No allocation in Update (RingBuffer + precomputed Hanning weights)
/// - True O(N log N) radix-2 FFT with bit-reversal permutation
/// - Pre-allocated work arrays for zero-allocation streaming
/// - Parabolic interpolation on peak bin for sub-bin accuracy
/// </remarks>
[SkipLocalsInit]
@@ -31,8 +35,10 @@ public sealed class Fft : AbstractBase
private readonly int _maxPeriod;
private readonly int _minBin;
private readonly int _maxBin;
private readonly double _twoPiOverN;
private readonly double[] _hanning;
private readonly int[] _bitRev;
private readonly double[] _workRe;
private readonly double[] _workIm;
private readonly RingBuffer _buffer;
[StructLayout(LayoutKind.Auto)]
@@ -44,7 +50,7 @@ public sealed class Fft : AbstractBase
/// <summary>
/// Initializes a new Fft indicator.
/// </summary>
/// <param name="windowSize">DFT window size in bars. Must be 32, 64, or 128. Default 64.</param>
/// <param name="windowSize">FFT window size in bars. Must be 32, 64, or 128. Default 64.</param>
/// <param name="minPeriod">Minimum detectable cycle period. Must be >= 2. Default 4.</param>
/// <param name="maxPeriod">Maximum detectable cycle period. Must be &lt;= windowSize/2. Default 32.</param>
public Fft(int windowSize = 64, int minPeriod = 4, int maxPeriod = 32)
@@ -67,19 +73,31 @@ public sealed class Fft : AbstractBase
_windowSize = windowSize;
_minPeriod = minPeriod;
_maxPeriod = maxPeriod;
_twoPiOverN = 2.0 * Math.PI / windowSize;
int log2N = Log2(windowSize);
// bin k corresponds to period N/k; k=minBin → period=N/minBin=maxPeriod, k=maxBin → period=N/maxBin=minPeriod
// bin k corresponds to period N/k
_minBin = Math.Max(1, windowSize / maxPeriod);
_maxBin = Math.Min(windowSize / 2, windowSize / minPeriod);
// Precompute Hanning window: w[n] = 0.5 - 0.5*cos(2π*n/N), n=0..N-1
// Precompute Hanning window: w[n] = 0.5 - 0.5*cos(2π*n/N)
double twoPiOverN = 2.0 * Math.PI / windowSize;
_hanning = new double[windowSize];
for (int n = 0; n < windowSize; n++)
{
_hanning[n] = 0.5 - 0.5 * Math.Cos(_twoPiOverN * n);
_hanning[n] = 0.5 - 0.5 * Math.Cos(twoPiOverN * n);
}
// Precompute bit-reversal permutation table
_bitRev = new int[windowSize];
for (int i = 0; i < windowSize; i++)
{
_bitRev[i] = BitReverse(i, log2N);
}
// Pre-allocate work arrays (zero allocation in hot path)
_workRe = new double[windowSize];
_workIm = new double[windowSize];
_buffer = new RingBuffer(windowSize);
Name = $"Fft({windowSize},{minPeriod},{maxPeriod})";
WarmupPeriod = windowSize;
@@ -90,10 +108,6 @@ public sealed class Fft : AbstractBase
/// <summary>
/// Initializes a new Fft indicator with source for event-based chaining.
/// </summary>
/// <param name="source">Source indicator for chaining</param>
/// <param name="windowSize">DFT window size. Must be 32, 64, or 128. Default 64.</param>
/// <param name="minPeriod">Minimum detectable period. Must be >= 2. Default 4.</param>
/// <param name="maxPeriod">Maximum detectable period. Must be &lt;= windowSize/2. Default 32.</param>
public Fft(ITValuePublisher source, int windowSize = 64, int minPeriod = 4, int maxPeriod = 32)
: this(windowSize, minPeriod, maxPeriod)
{
@@ -103,59 +117,155 @@ public sealed class Fft : AbstractBase
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void HandleUpdate(object? sender, in TValueEventArgs e) => Update(e.Value, e.IsNew);
/// <summary>
/// Computes floor(log2(n)) for powers of 2.
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static int Log2(int n)
{
int p = 0;
int x = n;
while (x > 1)
{
x >>= 1;
p++;
}
return p;
}
/// <summary>
/// Reverses the bits of x using 'bits' bit-width.
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static int BitReverse(int x, int bits)
{
int r = 0;
for (int i = 0; i < bits; i++)
{
r = (r << 1) | (x & 1);
x >>= 1;
}
return r;
}
/// <summary>
/// In-place iterative radix-2 Cooley-Tukey FFT.
/// </summary>
/// <param name="re">Real part array (modified in-place)</param>
/// <param name="im">Imaginary part array (modified in-place)</param>
/// <param name="n">Array length (must be power of 2)</param>
/// <param name="bitRev">Pre-computed bit-reversal table</param>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static void FftInPlace(double[] re, double[] im, int n, int[] bitRev)
{
// Bit-reversal permutation
for (int i = 0; i < n; i++)
{
int j = bitRev[i];
if (j > i)
{
(re[i], re[j]) = (re[j], re[i]);
(im[i], im[j]) = (im[j], im[i]);
}
}
// Cooley-Tukey butterfly stages
int len = 2;
while (len <= n)
{
int half = len >> 1;
double angStep = -2.0 * Math.PI / len;
for (int start = 0; start < n; start += len)
{
for (int k = 0; k < half; k++)
{
double angle = angStep * k;
double wr = Math.Cos(angle);
double wi = Math.Sin(angle);
int i0 = start + k;
int i1 = i0 + half;
double ur = re[i0];
double ui = im[i0];
double vr = re[i1];
double vi = im[i1];
// Twiddle: t = w * v
double tr = Math.FusedMultiplyAdd(vr, wr, -(vi * wi));
double ti = Math.FusedMultiplyAdd(vr, wi, vi * wr);
re[i0] = ur + tr;
im[i0] = ui + ti;
re[i1] = ur - tr;
im[i1] = ui - ti;
}
}
len <<= 1;
}
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private double ComputeDominantPeriod()
{
var span = _buffer.GetSpan();
int n = _windowSize;
double maxMag = 0.0;
int peakBin = _minBin;
double magBefore = 0.0;
double magAtPeak = 0.0;
double magAfter = 0.0;
// Fill work arrays: windowed data (oldest→newest), imag=0
for (int i = 0; i < n; i++)
{
_workRe[i] = span[i] * _hanning[i];
_workIm[i] = 0.0;
}
// Radix-2 FFT in-place
FftInPlace(_workRe, _workIm, n, _bitRev);
// Find peak magnitude in [minBin..maxBin]
double bestMag = -1.0;
int bestK = _minBin;
for (int k = _minBin; k <= _maxBin; k++)
{
double omegaK = _twoPiOverN * k;
double re = 0.0;
double im = 0.0;
for (int idx = 0; idx < n; idx++)
double mag = Math.FusedMultiplyAdd(_workRe[k], _workRe[k], _workIm[k] * _workIm[k]);
if (mag > bestMag)
{
// span[0]=oldest, span[n-1]=newest
// n=0 in DFT = current (newest): map DFT-n to span index (n-1-dftN)
// span[n-1-dftN]: dftN=0 → span[n-1] (newest), dftN=n-1 → span[0] (oldest)
double val = span[n - 1 - idx];
double xw = val * _hanning[idx];
double angle = omegaK * idx;
double cosA = Math.Cos(angle);
double sinA = Math.Sin(angle);
re = Math.FusedMultiplyAdd(xw, cosA, re);
im = Math.FusedMultiplyAdd(xw, -sinA, im);
}
double mag = Math.FusedMultiplyAdd(re, re, im * im);
if (mag > maxMag)
{
magBefore = magAtPeak;
magAfter = 0.0;
maxMag = mag;
magAtPeak = mag;
peakBin = k;
}
else if (peakBin > 0 && magAfter == 0.0)
{
magAfter = mag;
bestMag = mag;
bestK = k;
}
}
// Parabolic interpolation for sub-bin refinement
double denom = magBefore + 2.0 * maxMag + magAfter;
double shift = (denom > 0.0) ? (magBefore - magAfter) / denom : 0.0;
double dominantPeriod = (double)_windowSize / (peakBin + shift);
// Neighbor magnitudes for parabolic interpolation
double a, b, c;
b = bestMag;
if (bestK > _minBin)
{
a = Math.FusedMultiplyAdd(_workRe[bestK - 1], _workRe[bestK - 1],
_workIm[bestK - 1] * _workIm[bestK - 1]);
}
else
{
a = b;
}
if (bestK < _maxBin)
{
c = Math.FusedMultiplyAdd(_workRe[bestK + 1], _workRe[bestK + 1],
_workIm[bestK + 1] * _workIm[bestK + 1]);
}
else
{
c = b;
}
// Parabolic interpolation: shift = 0.5*(a-c)/(a - 2b + c)
double denom = a - 2.0 * b + c;
double shift = Math.Abs(denom) > 0.0 ? 0.5 * (a - c) / denom : 0.0;
double dominantPeriod = (double)_windowSize / (bestK + shift);
// Clamp to [minPeriod, maxPeriod]
return Math.Clamp(dominantPeriod, _minPeriod, _maxPeriod);
}
@@ -215,11 +325,6 @@ public sealed class Fft : AbstractBase
/// <summary>
/// Primes the indicator with historical values.
/// </summary>
/// <remarks>
/// Synthetic timestamps are generated by subtracting <c>step × source.Length</c>
/// from <see cref="DateTime.UtcNow"/>. For deterministic or replay-safe pipelines
/// use <see cref="Update(TValue, bool)"/> directly with explicit timestamps.
/// </remarks>
public override void Prime(ReadOnlySpan<double> source, TimeSpan? step = null)
{
TimeSpan interval = step ?? TimeSpan.FromSeconds(1);
@@ -239,8 +344,8 @@ public sealed class Fft : AbstractBase
}
/// <summary>
/// Computes dominant cycle period over a span of values using a sliding Hanning-windowed DFT.
/// Uses stackalloc for Hanning weights when windowSize &lt;= 64, otherwise ArrayPool.
/// Computes dominant cycle period over a span using sliding Hanning-windowed radix-2 FFT.
/// Uses stackalloc for work arrays when windowSize &lt;= 64, otherwise ArrayPool.
/// </summary>
public static void Batch(
ReadOnlySpan<double> src, Span<double> output,
@@ -271,15 +376,24 @@ public sealed class Fft : AbstractBase
throw new ArgumentException($"maxPeriod must be <= windowSize/2", nameof(maxPeriod));
}
int log2N = Log2(windowSize);
double twoPiOverN = 2.0 * Math.PI / windowSize;
int minBin = Math.Max(1, windowSize / maxPeriod);
int maxBin = Math.Min(windowSize / 2, windowSize / minPeriod);
double defaultPeriod = (minPeriod + maxPeriod) * 0.5;
double lastValid = defaultPeriod;
// Precompute Hanning window and bit-reversal table
const int StackallocThreshold = 64;
double[]? rentedW = null;
double[]? rentedH = null;
double[]? rentedRe = null;
double[]? rentedIm = null;
int[]? rentedBr = null;
scoped Span<double> hanning;
double[] workRe;
double[] workIm;
int[] bitRev;
if (windowSize <= StackallocThreshold)
{
@@ -287,15 +401,24 @@ public sealed class Fft : AbstractBase
}
else
{
rentedW = ArrayPool<double>.Shared.Rent(windowSize);
hanning = rentedW.AsSpan(0, windowSize);
rentedH = ArrayPool<double>.Shared.Rent(windowSize);
hanning = rentedH.AsSpan(0, windowSize);
}
// FFT work arrays (must be double[] for FftInPlace)
rentedRe = ArrayPool<double>.Shared.Rent(windowSize);
rentedIm = ArrayPool<double>.Shared.Rent(windowSize);
rentedBr = ArrayPool<int>.Shared.Rent(windowSize);
workRe = rentedRe;
workIm = rentedIm;
bitRev = rentedBr;
try
{
for (int n = 0; n < windowSize; n++)
{
hanning[n] = 0.5 - 0.5 * Math.Cos(twoPiOverN * n);
bitRev[n] = BitReverse(n, log2N);
}
for (int i = 0; i < src.Length; i++)
@@ -313,52 +436,49 @@ public sealed class Fft : AbstractBase
continue;
}
double maxMag = 0.0;
int peakBin = minBin;
double magBefore = 0.0;
double magAtPeak = 0.0;
double magAfter = 0.0;
// Fill work arrays with windowed data
for (int n = 0; n < windowSize; n++)
{
double v = src[i - windowSize + 1 + n];
if (!double.IsFinite(v))
{
v = lastValid;
}
workRe[n] = v * hanning[n];
workIm[n] = 0.0;
}
// Radix-2 FFT
FftInPlace(workRe, workIm, windowSize, bitRev);
// Find peak magnitude
double bestMag = -1.0;
int bestK = minBin;
for (int k = minBin; k <= maxBin; k++)
{
double omegaK = twoPiOverN * k;
double re = 0.0;
double im = 0.0;
for (int dftN = 0; dftN < windowSize; dftN++)
double mag = Math.FusedMultiplyAdd(workRe[k], workRe[k], workIm[k] * workIm[k]);
if (mag > bestMag)
{
// dftN=0 → newest (src[i]), dftN=windowSize-1 → oldest (src[start])
double v = src[i - dftN];
if (!double.IsFinite(v))
{
v = lastValid;
}
double xw = v * hanning[dftN];
double angle = omegaK * dftN;
re = Math.FusedMultiplyAdd(xw, Math.Cos(angle), re);
im = Math.FusedMultiplyAdd(xw, -Math.Sin(angle), im);
}
double mag = Math.FusedMultiplyAdd(re, re, im * im);
if (mag > maxMag)
{
magBefore = magAtPeak;
magAfter = 0.0;
maxMag = mag;
magAtPeak = mag;
peakBin = k;
}
else if (peakBin > 0 && magAfter == 0.0)
{
magAfter = mag;
bestMag = mag;
bestK = k;
}
}
double denom = magBefore + 2.0 * maxMag + magAfter;
double shift = (denom > 0.0) ? (magBefore - magAfter) / denom : 0.0;
double dominant = (double)windowSize / (peakBin + shift);
// Neighbor magnitudes for parabolic interpolation
double a = bestK > minBin
? Math.FusedMultiplyAdd(workRe[bestK - 1], workRe[bestK - 1],
workIm[bestK - 1] * workIm[bestK - 1])
: bestMag;
double c = bestK < maxBin
? Math.FusedMultiplyAdd(workRe[bestK + 1], workRe[bestK + 1],
workIm[bestK + 1] * workIm[bestK + 1])
: bestMag;
double denom = a - 2.0 * bestMag + c;
double shift = Math.Abs(denom) > 0.0 ? 0.5 * (a - c) / denom : 0.0;
double dominant = (double)windowSize / (bestK + shift);
double clamped = Math.Clamp(dominant, minPeriod, maxPeriod);
lastValid = clamped;
output[i] = clamped;
@@ -366,10 +486,13 @@ public sealed class Fft : AbstractBase
}
finally
{
if (rentedW != null)
if (rentedH != null)
{
ArrayPool<double>.Shared.Return(rentedW);
ArrayPool<double>.Shared.Return(rentedH);
}
ArrayPool<double>.Shared.Return(rentedRe);
ArrayPool<double>.Shared.Return(rentedIm);
ArrayPool<int>.Shared.Return(rentedBr);
}
}
lib/numerics/fft/Fft.md
+71 -52
View File
@@ -6,46 +6,51 @@
| **Inputs** | Source (close) |
| **Parameters** | `windowSize` (default 64), `minPeriod` (default 4), `maxPeriod` (default 32) |
| **Outputs** | Single series (Fft) |
| **Output range** | Varies (see docs) |
| **Warmup** | 1 bar |
| **Output range** | [minPeriod, maxPeriod] |
| **Warmup** | windowSize bars |
### TL;DR
- The FFT indicator computes the dominant cycle period in a price series using a Discrete Fourier Transform with a Hanning window.
- Parameterized by `windowsize` (default 64), `minperiod` (default 4), `maxperiod` (default 32).
- Output range: Varies (see docs).
- Requires 1 bar of warmup before first valid output (IsHot = true).
- Validated against TA-Lib, Skender, and Tulip reference implementations where available.
- The FFT indicator computes the dominant cycle period in a price series using a radix-2 Cooley-Tukey Fast Fourier Transform with a Hanning window.
- Parameterized by `windowSize` (default 64), `minPeriod` (default 4), `maxPeriod` (default 32).
- Output range: [minPeriod, maxPeriod] bars.
- Requires windowSize bars of warmup before first valid output (IsHot = true).
- True $O(N \log N)$ radix-2 FFT with bit-reversal permutation and Cooley-Tukey butterflies.
The FFT indicator computes the dominant cycle period in a price series using a Discrete Fourier Transform with a Hanning window. Rather than outputting frequency-domain magnitudes, it returns the estimated dominant cycle period in bars, making it directly usable as an adaptive period input for other indicators. The implementation uses a brute-force DFT over a constrained frequency band (not a radix-2 FFT), with parabolic interpolation on the magnitude spectrum to achieve sub-bin frequency resolution. With window sizes of 32, 64, or 128 and $O(N \cdot N/2)$ complexity per bar, the indicator trades computational cost for precise cycle detection within user-specified period bounds.
The FFT indicator computes the dominant cycle period in a price series using a true radix-2 Cooley-Tukey Fast Fourier Transform with a Hanning window. Rather than outputting frequency-domain magnitudes, it returns the estimated dominant cycle period in bars, making it directly usable as an adaptive period input for other indicators. The implementation uses an in-place iterative radix-2 FFT with bit-reversal permutation and Cooley-Tukey butterfly operations, achieving $O(N \log N)$ complexity. Parabolic interpolation on the magnitude spectrum provides sub-bin frequency resolution. With window sizes restricted to powers of two (32, 64, or 128), the indicator achieves precise cycle detection within user-specified period bounds with pre-allocated work arrays for zero-allocation streaming.
## Historical Context
The Fourier transform, formalized by Joseph Fourier (1822), decomposes any periodic signal into sinusoidal components. The Fast Fourier Transform algorithm (Cooley and Tukey, 1965) reduced the DFT from $O(N^2)$ to $O(N \log N)$, enabling real-time spectral analysis. However, for the small window sizes used in financial cycle detection (32-128 samples), the asymptotic advantage of FFT over DFT is minimal, and the DFT avoids the power-of-two length constraint.
The Fourier transform, formalized by Joseph Fourier (1822), decomposes any periodic signal into sinusoidal components. The Fast Fourier Transform algorithm, published by James Cooley and John Tukey in 1965, reduced the DFT from $O(N^2)$ to $O(N \log N)$ by recursively decomposing the DFT into smaller sub-problems using the "butterfly" operation pattern. The radix-2 variant requires power-of-two input lengths and uses bit-reversal permutation followed by iterative butterfly stages.
John Ehlers pioneered the application of spectral analysis to financial markets in the 1990s and 2000s, using DFT-based cycle measurement to create adaptive indicators. His work demonstrated that financial time series contain quasi-periodic cycles with time-varying periods, typically in the 6-40 bar range. The dominant cycle period, extracted via spectral peak detection, can drive adaptive moving averages (MAMA, FAMA), adaptive RSI, and other indicators that benefit from knowing the current market rhythm.
John Ehlers pioneered the application of spectral analysis to financial markets in the 1990s and 2000s, using FFT-based cycle measurement to create adaptive indicators. His work demonstrated that financial time series contain quasi-periodic cycles with time-varying periods, typically in the 6-40 bar range. The dominant cycle period, extracted via spectral peak detection, can drive adaptive moving averages (MAMA, FAMA), adaptive RSI, and other indicators that benefit from knowing the current market rhythm.
The Hanning window (also called Hann window, after Julius von Hann) is applied to reduce spectral leakage. Without windowing, the sharp truncation of a finite data segment creates artificial high-frequency components that contaminate the spectrum. The Hanning window tapers the data to zero at both ends, suppressing sidelobes at the cost of slightly wider main lobes (reduced frequency resolution).
## Architecture and Physics
The computation pipeline has four stages:
The computation pipeline has five stages:
**Stage 1: Windowed DFT** computes the real and imaginary components of the Fourier coefficients for frequency bins $k$ ranging from `minBin` to `maxBin`:
**Stage 1: Windowing** applies the Hanning window to the rolling price buffer:
$$X[k] = \sum_{n=0}^{N-1} x[n] \cdot w[n] \cdot e^{-j 2\pi k n / N}$$
$$x_w[n] = x[n] \cdot w[n], \quad w[n] = 0.5 - 0.5\cos\!\left(\frac{2\pi n}{N}\right)$$
where $w[n] = 0.5 - 0.5\cos(2\pi n/N)$ is the Hanning window. Only bins corresponding to periods in `[minPeriod, maxPeriod]` are evaluated, reducing computation.
**Stage 2: Bit-reversal permutation** reorders the windowed data according to the bit-reversed indices, preparing for in-place butterfly computation. The permutation table is pre-computed in the constructor.
**Stage 2: Power spectrum peak** finds the bin $k^*$ with maximum squared magnitude $|X[k]|^2 = \text{Re}^2 + \text{Im}^2$. During the search, the magnitudes of the bins adjacent to the peak (one before, one after) are captured for interpolation.
**Stage 3: Cooley-Tukey butterflies** perform $\log_2(N)$ stages of butterfly operations. Each stage $s$ processes pairs of elements separated by $2^{s-1}$ positions, combining them with twiddle factors:
**Stage 3: Parabolic interpolation** refines the peak location using a three-point parabola fit on the magnitudes at bins $k^*-1$, $k^*$, $k^*+1$:
$$\begin{aligned}
X[i_0] &\leftarrow X[i_0] + W_N^k \cdot X[i_1] \\
X[i_1] &\leftarrow X[i_0] - W_N^k \cdot X[i_1]
\end{aligned}$$
$$\delta = \frac{M_{k^*-1} - M_{k^*+1}}{M_{k^*-1} + 2 M_{k^*} + M_{k^*+1}}$$
where $W_N^k = e^{-j2\pi k/N}$ is the twiddle factor.
The refined dominant period is $N / (k^* + \delta)$.
**Stage 4: Peak detection with parabolic interpolation** finds the bin $k^*$ with maximum squared magnitude in `[minBin, maxBin]`, then refines using a three-point parabolic fit:
**Stage 4: Clamping** ensures the output stays within `[minPeriod, maxPeriod]`.
$$\delta = \frac{0.5 \cdot (P[k^*-1] - P[k^*+1])}{P[k^*-1] - 2P[k^*] + P[k^*+1]}$$
**Stage 5: Period extraction and clamping** converts the refined bin index to period: $T = N / (k^* + \delta)$, clamped to `[minPeriod, maxPeriod]`.
**Window size trade-offs**: $N = 32$ gives coarse resolution (period bins spaced ~1 bar apart) but fast response; $N = 128$ gives fine resolution (~0.25 bar spacing) but sluggish adaptation. The default $N = 64$ balances resolution and responsiveness.
@@ -55,6 +60,12 @@ The **Discrete Fourier Transform** for $N$ samples:
$$X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j 2\pi k n / N}, \quad k = 0, 1, \ldots, N-1$$
**Radix-2 Cooley-Tukey decomposition** splits the DFT into even and odd indexed sub-problems:
$$X[k] = \sum_{r=0}^{N/2-1} x[2r] \cdot W_{N/2}^{kr} + W_N^k \sum_{r=0}^{N/2-1} x[2r+1] \cdot W_{N/2}^{kr}$$
This recursion, applied iteratively with bit-reversal permutation, achieves $O(N \log N)$ complexity.
**Hanning window**:
$$w[n] = 0.5 - 0.5\cos\!\left(\frac{2\pi n}{N}\right)$$
@@ -69,7 +80,7 @@ $$k_{\min} = \max\!\left(1,\; \left\lfloor\frac{N}{T_{\max}}\right\rfloor\right)
**Parabolic interpolation** for sub-bin precision:
$$\hat{k} = k^* + \frac{P[k^*-1] - P[k^*+1]}{P[k^*-1] + 2P[k^*] + P[k^*+1]}$$
$$\hat{k} = k^* + \frac{0.5 \cdot (P[k^*-1] - P[k^*+1])}{P[k^*-1] - 2P[k^*] + P[k^*+1]}$$
$$T_{\text{dominant}} = \frac{N}{\hat{k}}$$
@@ -78,28 +89,33 @@ $$T_{\text{dominant}} = \frac{N}{\hat{k}}$$
```
FFT(source, windowSize, minPeriod, maxPeriod):
N = windowSize
twoPiOverN = 2 * pi / N
minBin = max(1, N / maxPeriod)
maxBin = min(N/2, N / minPeriod)
// Stage 1: Apply Hanning window
for n = 0 to N-1:
workRe[n] = source[n] * hanning[n]
workIm[n] = 0
maxMag = 0; peakBin = 0
for k = minBin to maxBin:
re = 0; im = 0
for n = 0 to N-1:
w = 0.5 - 0.5 * cos(twoPiOverN * n) // Hanning
xw = source[n] * w
angle = twoPiOverN * k * n
re += xw * cos(angle)
im -= xw * sin(angle)
mag = re*re + im*im
if mag > maxMag:
track neighbor magnitudes
maxMag = mag; peakBin = k
// Stage 2: Bit-reversal permutation
for i = 0 to N-1:
j = bitReverse(i)
if j > i: swap(workRe[i], workRe[j])
// Parabolic interpolation
shift = (magBefore - magAfter) / (magBefore + 2*maxMag + magAfter)
dominantPeriod = N / (peakBin + shift)
return clamp(dominantPeriod, minPeriod, maxPeriod)
// Stage 3: Cooley-Tukey butterflies
len = 2
while len <= N:
half = len / 2
angStep = -2π / len
for start = 0 to N-1 step len:
for k = 0 to half-1:
w = exp(j * angStep * k)
butterfly(workRe, workIm, start+k, start+k+half, w)
len *= 2
// Stage 4: Peak detection + interpolation
peakBin = argmax |X[k]|² for k in [minBin..maxBin]
shift = 0.5*(P[k-1] - P[k+1]) / (P[k-1] - 2*P[k] + P[k+1])
// Stage 5: Period extraction
return clamp(N / (peakBin + shift), minPeriod, maxPeriod)
```
@@ -107,34 +123,37 @@ FFT(source, windowSize, minPeriod, maxPeriod):
### Operation Count (Streaming Mode)
FFT (DFT dominant cycle detector) evaluates B frequency bins, each requiring N multiply-accumulates — O(N*B) per bar.
FFT (radix-2 Cooley-Tukey) performs N/2 butterflies per stage across log₂(N) stages — O(N log N) per bar.
| Operation | Count | Cost (cycles) | Subtotal |
| :--- | :---: | :---: | :---: |
| Hanning window multiply | N | 2 cy | ~2N cy |
| DFT inner loop (B bins * N samples) | B*N | 4 cy | ~4*N*B cy |
| cos/sin evaluation (precomputed table) | 2*B*N | 0 cy | ~0 cy |
| Magnitude comparison + peak track | B | 2 cy | ~2B cy |
| Parabolic interpolation (3 points) | 1 | 5 cy | ~5 cy |
| **Total (N=64, B=10)** | **O(N*B)** | — | **~2617 cy** |
| Bit-reversal permutation | N | 1 cy | ~N cy |
| Butterfly operations (log₂N stages × N/2) | N/2 × log₂N | 8 cy | ~4N·log₂N cy |
| cos/sin per butterfly | N/2 × log₂N | 14 cy | ~7N·log₂N cy |
| Magnitude search (B bins) | B | 4 cy | ~4B cy |
| Parabolic interpolation | 1 | 10 cy | ~10 cy |
| **Total (N=64, B=10)** | **O(N log N)** | — | **~4362 cy** |
O(N*B) per bar where B = active frequency bins. Precomputed sin/cos tables eliminate transcendental cost. Suitable for 1-minute+ timeframes; not tick-data hot paths.
O(N log N) per bar. Pre-allocated work arrays ensure zero allocation in the hot path. Twiddle factor computation dominates; pre-computing sin/cos tables would reduce to ~2500 cy.
### Batch Mode (SIMD Analysis)
| Operation | Vectorizable? | Notes |
| :--- | :---: | :--- |
| Hanning window application | Yes | Vector multiply with precomputed weights |
| DFT inner dot product | Yes | FMA with sin/cos table lookup |
| Magnitude squared | Yes | Vector FMA (re^2 + im^2) |
| Bit-reversal permutation | No | Random access pattern; scalar only |
| Butterfly multiply-add | Yes | Complex FMA operations on paired elements |
| Magnitude squared | Yes | Vector FMA (re² + im²) |
| Peak search | Partial | Max reduction; SIMD-friendly |
Strong batch SIMD: inner dot products are FMA-vectorizable. AVX2 processes 4 complex outputs per 2 cycles. Expected 3-4× speedup for N=64.
Moderate SIMD potential: butterfly FMA operations are vectorizable within each stage. Bit-reversal permutation is inherently scalar. Expected 2× speedup over scalar for N=64.
## Resources
- Cooley, J.W. & Tukey, J.W. "An Algorithm for the Machine Calculation of Complex Fourier Series." Mathematics of Computation, 1965.
- Cooley, J.W. & Tukey, J.W. "An Algorithm for the Machine Calculation of Complex Fourier Series." *Mathematics of Computation*, 1965.
- Ehlers, J.F. "Cycle Analytics for Traders." Wiley, 2013.
- Ehlers, J.F. "Rocket Science for Traders." Wiley, 2001.
- Harris, F.J. "On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform." Proc. IEEE, 1978.
- Harris, F.J. "On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform." *Proc. IEEE*, 1978.
- Oppenheim, A.V. & Schafer, R.W. "Discrete-Time Signal Processing." 3rd edition, Pearson, 2010.
- PineScript reference: [`fft.pine`](fft.pine)