From 9582028f3b9de0e79a5c1ff2bc0e4f587747263f Mon Sep 17 00:00:00 2001 From: Justin Purnell Date: Tue, 7 Apr 2026 06:37:33 -1000 Subject: [PATCH 1/2] =?UTF-8?q?Fix=20AccelerateFFTBackend=204=C3=97=20powe?= =?UTF-8?q?r=20scaling=20bug?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit vDSP_fft_zripD returns FFT outputs scaled by 2 vs the textbook DFT formula (vDSP convention for packed real-input FFT). Squaring magnitudes therefore produced values 4× the textbook |X[k]|² on Darwin only, breaking absolute-power analysis (Parseval's theorem, PSD integration, band power, etc.). The pre-existing cross-backend test only checked peak bin location, not absolute magnitudes, so the discrepancy was invisible. The new PowerSpectralDensityTests suite (added in the next commit) is the lever that surfaced it. Fix: apply ×0.25 correction to all DC, Nyquist, and typical bins in AccelerateFFTBackend.powerSpectrum. Both backends now produce identical absolute power values to within machine precision. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../BusinessMath/Streaming/FFTBackend.swift | 20 +++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/Sources/BusinessMath/Streaming/FFTBackend.swift b/Sources/BusinessMath/Streaming/FFTBackend.swift index a214d2bd..ddb51477 100644 --- a/Sources/BusinessMath/Streaming/FFTBackend.swift +++ b/Sources/BusinessMath/Streaming/FFTBackend.swift @@ -246,22 +246,34 @@ public struct AccelerateFFTBackend: FFTBackend, Sendable { } } - // Compute power spectrum from results + // Compute power spectrum from results. + // + // **vDSP scaling correction:** `vDSP_fft_zripD` returns FFT outputs + // scaled by 2 vs the textbook DFT formula (vDSP convention for + // packed real-input FFT). Squaring magnitudes therefore yields + // values 4× the textbook |X[k]|². We divide by 4 to match the + // mathematical definition that `PureSwiftFFTBackend` produces, so + // both backends are interchangeable for absolute-power analysis + // (Parseval, PSD, band integration, etc.). + // + // Without this correction, downstream PSD computation would be + // off by 4× on Darwin only, breaking Parseval's theorem. let numBins = n / 2 + 1 var power = [Double](repeating: 0.0, count: numBins) + let vdspScalingCorrection = 0.25 // = 1/4 // DC component - power[0] = realPart[0] * realPart[0] + power[0] = realPart[0] * realPart[0] * vdspScalingCorrection // Nyquist component (packed in imagPart[0] for real FFT) if numBins > 1 { - power[numBins - 1] = imagPart[0] * imagPart[0] + power[numBins - 1] = imagPart[0] * imagPart[0] * vdspScalingCorrection } // Other bins for k in 1..<(n / 2) { if k < numBins { - power[k] = realPart[k] * realPart[k] + imagPart[k] * imagPart[k] + power[k] = (realPart[k] * realPart[k] + imagPart[k] * imagPart[k]) * vdspScalingCorrection } } From 1abbc70f7f286526866d5adc7e8293f494d8111a Mon Sep 17 00:00:00 2001 From: Justin Purnell Date: Tue, 7 Apr 2026 06:38:33 -1000 Subject: [PATCH 2/2] Add powerSpectralDensity to FFTBackend (v2.1.1) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds a normalized one-sided power spectral density method to the FFTBackend protocol so downstream consumers (BioFeedbackKit, etc.) get physically meaningful spectral values directly without reinventing normalization. API additions: - FFTBackend.powerSpectralDensity(_:sampleRate:) — one-sided PSD in units²/Hz. Default implementation in an extension; all backends inherit for free. - FFTBackend.powerSpectralDensityBins(_:sampleRate:) — convenience that pairs each PSD value with its center frequency in Hz. - PSDBin value type — Sendable, Equatable, (frequency: Double, power: Double). Normalization correctness: - One-sided spectrum: DC and Nyquist bins use edge factor 1/(M·fs); typical bins use 2/(M·fs). - Uses the UNPADDED signal length M, not the internally zero-padded length N. This ensures the PSD integral equals the time-domain variance per Parseval's theorem regardless of input length. Previously every downstream consumer had to know about this gotcha. Testing: - 12 new tests in PowerSpectralDensityTests.swift covering Parseval on pure sines, two-tone signals, the M=50→64 zero-padding edge case, DC and Nyquist bin handling, cross-backend equivalence (PureSwift vs Accelerate, both now satisfy Parseval after the previous commit's scaling fix), empty/invalid input handling, and PSDBin convenience. - Validation playground at Tests/Validation/PSD_Validation.swift — standalone hand-rolled implementation that verifies the formulas against Parseval's theorem independent of BusinessMath, producing the exact values the test suite asserts. - Full BusinessMath suite: 4720/4720 passing, zero compiler warnings. Purely additive — no breaking changes. Existing powerSpectrum(_:) signature is unchanged. Co-Authored-By: Claude Opus 4.6 (1M context) --- CHANGELOG.md | 51 ++++ .../BusinessMath/Streaming/FFTBackend.swift | 123 ++++++++ .../PowerSpectralDensityTests.swift | 272 ++++++++++++++++++ Tests/Validation/PSD_Validation.swift | 258 +++++++++++++++++ 4 files changed, 704 insertions(+) create mode 100644 Tests/BusinessMathTests/StreamingTests/PowerSpectralDensityTests.swift create mode 100644 Tests/Validation/PSD_Validation.swift diff --git a/CHANGELOG.md b/CHANGELOG.md index b96f91e0..249e950f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,57 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## BusinessMath Library +### [2.1.1] - 2026-04-06 + +**Version 2.1.1** adds normalized power spectral density (PSD) to the FFT layer +and fixes a pre-existing 4× scaling bug in `AccelerateFFTBackend.powerSpectrum` +that was uncovered while implementing PSD. Driven by the BioFeedbackKit project, +which needs physically meaningful spectral magnitudes (ms² for HRV LF/HF +analysis) without every consumer reinventing FFT normalization. + +#### Added + +- **`FFTBackend.powerSpectralDensity(_:sampleRate:)`** — new protocol method + returning a one-sided PSD in `units²/Hz`. The integral over frequency + equals the time-domain variance of a zero-mean input (Parseval's theorem), + so band powers come out in physical units directly. Default implementation + in an extension; all conformers (`PureSwiftFFTBackend`, `AccelerateFFTBackend`) + inherit it for free. +- **Normalization correctness:** the PSD method uses the **unpadded** signal + length `M` for normalization, not the internally zero-padded length `N`, + so PSD integrals remain physically meaningful regardless of input length. + Previously, every downstream consumer needed to know about this gotcha; + now BusinessMath handles it. +- **`PSDBin` value type** — pairs each PSD value with its center frequency in + Hz, returned by the convenience method `powerSpectralDensityBins(_:sampleRate:)`. +- **12 new tests** in `Tests/BusinessMathTests/StreamingTests/PowerSpectralDensityTests.swift` + covering Parseval's theorem on multiple fixtures, the M-vs-N normalization + edge case (M=50 padded to 64), DC and Nyquist edge factors, cross-backend + equivalence, and the `PSDBin` convenience. +- **Validation playground** at `Tests/Validation/PSD_Validation.swift` — + standalone hand-rolled implementation that verifies the PSD formulas + against Parseval's theorem before any package code runs. + +#### Fixed + +- **`AccelerateFFTBackend.powerSpectrum` 4× scaling bug.** `vDSP_fft_zripD` + returns FFT outputs scaled by 2 vs the textbook DFT formula (vDSP + convention for packed real-input FFT). Squaring magnitudes therefore + produced values 4× the textbook `|X[k]|²` on Darwin only. The existing + cross-backend test (`Accelerate FFT: matches Pure Swift within tolerance`) + only checked peak bin location, not absolute values, so the discrepancy + was invisible to it. Fix: apply a `× 0.25` correction in the power + computation. Both backends now satisfy Parseval's theorem and produce + identical PSD values to within `1e-9` relative tolerance. + +#### Notes + +- This release is **purely additive at the public API level**. The existing + `powerSpectrum(_:)` method signature is unchanged. Consumers that were + comparing absolute spectrum magnitudes from `AccelerateFFTBackend` against + external references will see corrected (smaller by 4×) values — but those + comparisons were previously wrong, and any such consumers should update. + ### [2.1.0] - 2026-04-06 **Version 2.1.0** is a stability and quality release that fixes a CI crash (SIGABRT), eliminates all compiler warnings, hardens concurrency safety, and adds Thread Sanitizer to the CI pipeline. diff --git a/Sources/BusinessMath/Streaming/FFTBackend.swift b/Sources/BusinessMath/Streaming/FFTBackend.swift index ddb51477..3f4c3d80 100644 --- a/Sources/BusinessMath/Streaming/FFTBackend.swift +++ b/Sources/BusinessMath/Streaming/FFTBackend.swift @@ -45,6 +45,129 @@ public protocol FFTBackend: Sendable { /// - Returns: Power spectrum values of length `N/2 + 1` where N is the /// zero-padded signal length. func powerSpectrum(_ signal: [Double]) -> [Double] + + /// Compute the one-sided power spectral density (PSD) of a real-valued signal. + /// + /// The integral of the returned PSD over frequency equals the time-domain + /// variance of the input signal (Parseval's theorem). Values are in + /// `units²/Hz` where `units` is the unit of the input signal — for example, + /// HRV RR-interval data in milliseconds yields PSD in `ms²/Hz` and + /// integrated band power in `ms²`. + /// + /// **Normalization conventions:** + /// - One-sided spectrum: bins `1.. [Double] +} + +// MARK: - PSD Default Implementation + +extension FFTBackend { + + /// Default one-sided PSD implementation, derived from `powerSpectrum(_:)`. + /// + /// All conforming backends inherit this for free. Backends MAY override + /// for performance, but the default is correct and stable. + public func powerSpectralDensity( + _ signal: [Double], + sampleRate: Double + ) -> [Double] { + guard signal.isEmpty == false, sampleRate > 0 else { return [] } + + let unpaddedLength = signal.count + let raw = powerSpectrum(signal) + guard raw.isEmpty == false else { return [] } + + let nyquistBin = raw.count - 1 + // Typical bins: factor of 2 for the one-sided convention + let typicalFactor = 2.0 / (Double(unpaddedLength) * sampleRate) + // DC bin and Nyquist bin: no factor of 2 + let edgeFactor = 1.0 / (Double(unpaddedLength) * sampleRate) + + var psd = [Double](repeating: 0.0, count: raw.count) + psd[0] = raw[0] * edgeFactor + if nyquistBin > 0 { + psd[nyquistBin] = raw[nyquistBin] * edgeFactor + } + if nyquistBin > 1 { + for k in 1.. [PSDBin] { + let psd = powerSpectralDensity(signal, sampleRate: sampleRate) + guard psd.isEmpty == false else { return [] } + let paddedLength = (psd.count - 1) * 2 + guard paddedLength > 0 else { return [] } + let deltaF = sampleRate / Double(paddedLength) + return psd.enumerated().map { idx, value in + PSDBin(frequency: Double(idx) * deltaF, power: value) + } + } +} + +// MARK: - PSDBin + +/// A single power spectral density value paired with its center frequency. +/// +/// Returned by ``FFTBackend/powerSpectralDensityBins(_:sampleRate:)``. +public struct PSDBin: Sendable, Equatable { + + /// Center frequency of the bin, in Hz. + public let frequency: Double + + /// Power spectral density at this bin, in `units²/Hz`. + public let power: Double + + /// Creates a PSD bin with the given frequency and power. + public init(frequency: Double, power: Double) { + self.frequency = frequency + self.power = power + } } // MARK: - Pure Swift FFT Backend diff --git a/Tests/BusinessMathTests/StreamingTests/PowerSpectralDensityTests.swift b/Tests/BusinessMathTests/StreamingTests/PowerSpectralDensityTests.swift new file mode 100644 index 00000000..f423a001 --- /dev/null +++ b/Tests/BusinessMathTests/StreamingTests/PowerSpectralDensityTests.swift @@ -0,0 +1,272 @@ +// +// PowerSpectralDensityTests.swift +// BusinessMath +// +// Created by Justin Purnell on 2026-04-06. +// +// Tests for the v2.1.1 powerSpectralDensity additions to FFTBackend. +// Reference truth: Parseval's theorem — the integral of a one-sided PSD +// over frequency must equal the time-domain variance of the (zero-mean) +// input signal. +// + +import Testing +import Foundation +@testable import BusinessMath + +@Suite("Power Spectral Density (v2.1.1)") +struct PowerSpectralDensityTests { + + // MARK: - Helpers + + /// Subtract the mean from a signal so Parseval's theorem holds without + /// the DC term skewing comparisons. + static func meanRemoved(_ x: [Double]) -> [Double] { + let m = x.reduce(0, +) / Double(x.count) + return x.map { $0 - m } + } + + /// Sample variance with `n` denominator (matches Σx²/n for zero-mean signals). + static func populationVariance(_ x: [Double]) -> Double { + let m = x.reduce(0, +) / Double(x.count) + let sq = x.map { ($0 - m) * ($0 - m) } + return sq.reduce(0, +) / Double(x.count) + } + + /// Integrate a one-sided PSD over frequency: Σ PSD[k] · Δf where Δf = fs/N_padded. + static func parsevalIntegral(psd: [Double], sampleRate: Double) -> Double { + guard !psd.isEmpty else { return 0 } + let N = (psd.count - 1) * 2 + guard N > 0 else { return 0 } + let deltaF = sampleRate / Double(N) + return psd.reduce(0, +) * deltaF + } + + // MARK: - Parseval — pure sine wave (no padding) + + @Test("Parseval: pure sine, M=64 (power of 2, no padding)") + func parsevalPureSineNoPadding() { + let backend = PureSwiftFFTBackend() + let M = 64 + let fs = 64.0 + let f0 = 4.0 + let A = 1.0 + let signal = Self.meanRemoved((0.. otherMax) + } + + // MARK: - Smoke tests + + @Test("Empty signal returns empty PSD") + func emptySignalReturnsEmpty() { + let backend = PureSwiftFFTBackend() + #expect(backend.powerSpectralDensity([], sampleRate: 100.0).isEmpty) + } + + @Test("Non-positive sample rate returns empty PSD") + func nonPositiveSampleRateReturnsEmpty() { + let backend = PureSwiftFFTBackend() + let signal = [Double](repeating: 1.0, count: 16) + #expect(backend.powerSpectralDensity(signal, sampleRate: 0.0).isEmpty) + #expect(backend.powerSpectralDensity(signal, sampleRate: -1.0).isEmpty) + } + + // MARK: - Backend equivalence (Darwin only) + + #if canImport(Accelerate) + @Test("Parseval holds for PureSwiftFFTBackend on a two-tone signal") + func parsevalPureSwift() { + let backend = PureSwiftFFTBackend() + let M = 256 + let fs = 256.0 + let signal = Self.meanRemoved((0.. Int { + guard n > 0 else { return 1 } + var p = 1 + while p < n { p *= 2 } + return p +} + +func bitReverse(_ x: Int, bits: Int) -> Int { + var r = 0 + var v = x + for _ in 0..>= 1 + } + return r +} + +/// In-place radix-2 FFT. `real`/`imag` arrays must have length n, a power of 2. +func fft(_ real: inout [Double], _ imag: inout [Double], _ n: Int) { + let logN = Int(log2(Double(n))) + for i in 0.. [Double] { + guard !signal.isEmpty else { return [] } + let n = nextPowerOf2(signal.count) + var real = [Double](repeating: 0, count: n) + var imag = [Double](repeating: 0, count: n) + for i in 0.. [Double] { + guard !signal.isEmpty, sampleRate > 0 else { return [] } + let M = signal.count + let raw = rawPowerSpectrum(signal) + guard !raw.isEmpty else { return [] } + let nyquistBin = raw.count - 1 + let typicalFactor = 2.0 / (Double(M) * sampleRate) + let edgeFactor = 1.0 / (Double(M) * sampleRate) + var out = [Double](repeating: 0, count: raw.count) + out[0] = raw[0] * edgeFactor + if nyquistBin > 0 { + out[nyquistBin] = raw[nyquistBin] * edgeFactor + } + if nyquistBin > 1 { + for k in 1.. Double { + guard !x.isEmpty else { return 0 } + let mean = x.reduce(0, +) / Double(x.count) + let sq = x.map { ($0 - mean) * ($0 - mean) } + return sq.reduce(0, +) / Double(x.count) +} + +func parsevalIntegral(psd: [Double], sampleRate: Double) -> Double { + guard !psd.isEmpty else { return 0 } + let N = (psd.count - 1) * 2 + guard N > 0 else { return 0 } + let deltaF = sampleRate / Double(N) + return psd.reduce(0, +) * deltaF +} + +// MARK: - Test fixtures + +print("=== PSD Normalization Validation ===") +print() + +// MARK: Fixture 1 — Pure sine wave, M = 64 (power of 2, no padding) +do { + print("--- Fixture 1: Pure sine, M=64, fs=64 Hz, f0=4 Hz, A=1 ---") + let M = 64 + let fs = 64.0 + let f0 = 4.0 + let A = 1.0 + let x = (0.. [Double] { + let m = x.reduce(0, +) / Double(x.count) + return x.map { $0 - m } + } + let xFull = meanRemoved((0..<64).map { i in A * sin(2.0 * .pi * f0 * Double(i) / fs) }) + let xShort = meanRemoved(Array((0..<50).map { i in A * sin(2.0 * .pi * f0 * Double(i) / fs) })) + let varFull = variance(xFull) + let varShort = variance(xShort) + let pFull = psd(xFull, sampleRate: fs) + let pShort = psd(xShort, sampleRate: fs) + let intFull = parsevalIntegral(psd: pFull, sampleRate: fs) + let intShort = parsevalIntegral(psd: pShort, sampleRate: fs) + print("xFull (M=64): variance=\(varFull), PSD integral=\(intFull), rel.err=\(abs(intFull-varFull)/varFull)") + print("xShort (M=50): variance=\(varShort), PSD integral=\(intShort), rel.err=\(abs(intShort-varShort)/varShort)") + print("Both should now satisfy Parseval at machine epsilon — proves M-vs-N normalization is correct.") + print() +} + +// MARK: Fixture 5 — DC signal +do { + print("--- Fixture 5: DC signal, M=16, fs=16 Hz, value=5 ---") + let M = 16 + let fs = 16.0 + let x = [Double](repeating: 5.0, count: M) + let timeVar = variance(x) // 0 + let p = psd(x, sampleRate: fs) + print("Time-domain variance: \(timeVar) (should be 0)") + print("PSD bin 0 (DC): \(p[0])") + print("PSD bin 1: \(p[1]) (should be 0 within numerical noise)") + print("PSD Nyquist: \(p[p.count - 1])") + print("Sum of non-DC bins: \(p[1...].reduce(0, +))") + print("Note: DC bin is non-zero (signal has DC), but variance is zero,") + print("which is exactly what Parseval expects: variance is the AC content.") + print() +} + +// MARK: Fixture 6 — Nyquist-only signal +do { + print("--- Fixture 6: Nyquist signal, alternating +1/-1, M=16, fs=16 Hz ---") + let M = 16 + let fs = 16.0 + let x = (0..