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 a214d2bd..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..<N/2` carry a factor of 2; the DC bin
+    ///   `0` and the Nyquist bin `N/2` do **not**.
+    /// - The normalization uses the **unpadded** signal length `M`, not the
+    ///   internally zero-padded length `N`. This ensures the PSD integral
+    ///   equals the time-domain variance regardless of input length.
+    /// - No window function is applied. Callers that want a tapered window
+    ///   (Hann, Hamming, etc.) must apply it before calling, and compensate
+    ///   the result by dividing by the window's noise-equivalent bandwidth
+    ///   `(1/M) · Σ w[m]²`.
+    /// - For Parseval to hold exactly, the input should be zero-mean (apply
+    ///   mean removal before calling). Otherwise the DC bin reflects the
+    ///   signal's mean and the integral exceeds the variance by `mean²`.
+    ///
+    /// ## Example
+    /// ```swift
+    /// let backend = FFTBackendSelector.selectBackend()
+    /// let mean = signal.reduce(0, +) / Double(signal.count)
+    /// let zeroMean = signal.map { $0 - mean }
+    /// let psd = backend.powerSpectralDensity(zeroMean, sampleRate: 4.0)
+    /// ```
+    ///
+    /// - Parameters:
+    ///   - signal: Real-valued input signal. Should be mean-removed (and
+    ///     optionally windowed) before calling. Internally zero-padded to
+    ///     the next power of 2 for the FFT.
+    ///   - sampleRate: Sample rate in Hz. Must be positive.
+    /// - Returns: One-sided PSD bins of length `N/2 + 1` where `N` is the
+    ///   zero-padded length. Returns an empty array for an empty signal or
+    ///   non-positive sample rate.
+    func powerSpectralDensity(_ signal: [Double], sampleRate: Double) -> [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..<nyquistBin {
+                psd[k] = raw[k] * typicalFactor
+            }
+        }
+        return psd
+    }
+
+    /// One-sided PSD with each bin labeled by its center frequency in Hz.
+    ///
+    /// Convenience that pairs ``powerSpectralDensity(_:sampleRate:)`` values
+    /// with their bin frequencies, so downstream code doesn't need to
+    /// recompute `Δf = sampleRate / N_padded`.
+    ///
+    /// - Parameters:
+    ///   - signal: Real-valued input signal. See ``powerSpectralDensity(_:sampleRate:)``
+    ///     for preparation requirements.
+    ///   - sampleRate: Sample rate in Hz.
+    /// - Returns: An array of ``PSDBin`` values. Empty for empty input.
+    public func powerSpectralDensityBins(
+        _ signal: [Double],
+        sampleRate: Double
+    ) -> [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
@@ -246,22 +369,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
             }
         }
 
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..<M).map { i in
+            A * sin(2.0 * .pi * f0 * Double(i) / fs)
+        })
+
+        let timeVar = Self.populationVariance(signal)
+        let psd = backend.powerSpectralDensity(signal, sampleRate: fs)
+        let integral = Self.parsevalIntegral(psd: psd, sampleRate: fs)
+
+        #expect(abs(integral - timeVar) / timeVar < 1e-12)
+        // Sanity: also matches the analytic A²/2 for a pure sine
+        #expect(abs(timeVar - A * A / 2) < 1e-12)
+    }
+
+    // MARK: - Parseval — multiple tones
+
+    @Test("Parseval: two-tone signal, variances add")
+    func parsevalTwoTones() {
+        let backend = PureSwiftFFTBackend()
+        let M = 256
+        let fs = 256.0
+        let f1 = 8.0
+        let f2 = 32.0
+        let A1 = 2.0
+        let A2 = 3.0
+        let signal = Self.meanRemoved((0..<M).map { i in
+            let t = Double(i) / fs
+            return A1 * sin(2.0 * .pi * f1 * t) + A2 * sin(2.0 * .pi * f2 * t)
+        })
+
+        let timeVar = Self.populationVariance(signal)
+        let analyticVar = A1 * A1 / 2 + A2 * A2 / 2
+        let psd = backend.powerSpectralDensity(signal, sampleRate: fs)
+        let integral = Self.parsevalIntegral(psd: psd, sampleRate: fs)
+
+        #expect(abs(timeVar - analyticVar) < 1e-12)
+        #expect(abs(integral - timeVar) / timeVar < 1e-12)
+    }
+
+    // MARK: - The critical M-vs-N normalization test
+
+    @Test("Parseval with zero-padding: M=50 padded to 64 still satisfies Parseval")
+    func parsevalZeroPaddingEquivalence() {
+        let backend = PureSwiftFFTBackend()
+        let fs = 64.0
+        let f0 = 8.0
+        let A = 1.5
+
+        // M=50 — NOT a power of 2, will be padded internally to 64
+        let signal = Self.meanRemoved((0..<50).map { i in
+            A * sin(2.0 * .pi * f0 * Double(i) / fs)
+        })
+
+        let timeVar = Self.populationVariance(signal)
+        let psd = backend.powerSpectralDensity(signal, sampleRate: fs)
+        let integral = Self.parsevalIntegral(psd: psd, sampleRate: fs)
+
+        // This is the test that proves the M-vs-N fix:
+        // if normalization used the padded length N=64 instead of M=50,
+        // PSD values would be too small by a factor of 50/64 and the
+        // integral would be ~22% off, not at machine epsilon.
+        #expect(abs(integral - timeVar) / timeVar < 1e-12)
+    }
+
+    // MARK: - Output bin count
+
+    @Test("PSD output length matches N_padded/2 + 1")
+    func outputLengthMatchesPaddedNyquist() {
+        let backend = PureSwiftFFTBackend()
+        // M=100 → padded to 128 → 65 bins
+        let signal = [Double](repeating: 0.0, count: 100)
+        let psd = backend.powerSpectralDensity(signal, sampleRate: 100.0)
+        #expect(psd.count == 128 / 2 + 1)
+    }
+
+    // MARK: - DC bin (not doubled in one-sided convention)
+
+    @Test("DC bin uses edge factor (not doubled)")
+    func dcBinEdgeFactor() {
+        let backend = PureSwiftFFTBackend()
+        // Pure DC signal — all power should land in bin 0
+        let M = 16
+        let fs = 16.0
+        let signal = [Double](repeating: 5.0, count: M)
+        let psd = backend.powerSpectralDensity(signal, sampleRate: fs)
+
+        // Raw |X[0]|² for a constant signal of value c is (c·M)² = c²·M²
+        // Edge factor for DC: 1/(M·fs)
+        // So PSD[0] should equal c² · M² · 1/(M·fs) = c² · M / fs = 25 · 16 / 16 = 25
+        #expect(abs(psd[0] - 25.0) < 1e-12)
+
+        // All non-DC bins should be ~0
+        for k in 1..<psd.count {
+            #expect(abs(psd[k]) < 1e-12)
+        }
+    }
+
+    // MARK: - Nyquist bin (also not doubled)
+
+    @Test("Nyquist bin uses edge factor (not doubled)")
+    func nyquistBinEdgeFactor() {
+        let backend = PureSwiftFFTBackend()
+        // Alternating +1, -1 — all energy at exactly the Nyquist frequency
+        let M = 16
+        let fs = 16.0
+        let signal = (0..<M).map { i in i.isMultiple(of: 2) ? 1.0 : -1.0 }
+        let timeVar = Self.populationVariance(signal)  // = 1.0
+
+        let psd = backend.powerSpectralDensity(signal, sampleRate: fs)
+        let integral = Self.parsevalIntegral(psd: psd, sampleRate: fs)
+
+        // Variance is 1.0 for ±1 alternation
+        #expect(abs(timeVar - 1.0) < 1e-12)
+        // Parseval must hold
+        #expect(abs(integral - timeVar) < 1e-12)
+        // The Nyquist bin (last bin) should carry the power, not other bins
+        let nyq = psd.last ?? 0
+        let otherMax = psd.dropLast().max() ?? 0
+        #expect(nyq > 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..<M).map { i in
+            let t = Double(i) / fs
+            return 2.0 * sin(2.0 * .pi * 8.0 * t) + 3.0 * sin(2.0 * .pi * 32.0 * t)
+        })
+        let timeVar = Self.populationVariance(signal)
+        let psd = backend.powerSpectralDensity(signal, sampleRate: fs)
+        let integral = Self.parsevalIntegral(psd: psd, sampleRate: fs)
+        #expect(abs(integral - timeVar) / timeVar < 1e-12)
+    }
+
+    /// Cross-backend equivalence: PureSwift and Accelerate must produce
+    /// PSDs that agree to within numerical tolerance on the same input.
+    ///
+    /// This test was the lever that surfaced a pre-existing 4× scaling bug
+    /// in `AccelerateFFTBackend.powerSpectrum`: `vDSP_fft_zripD` returns
+    /// FFT outputs scaled by 2 vs the textbook DFT formula, and squaring
+    /// produced values 4× the correct `|X[k]|²`. The fix (a `× 0.25`
+    /// correction in the power computation) is part of v2.1.1 alongside
+    /// the new PSD method.
+    @Test("Cross-backend equivalence: PureSwift and Accelerate produce identical PSDs")
+    func crossBackendEquivalence() {
+        let pureSwift = PureSwiftFFTBackend()
+        let accelerate = AccelerateFFTBackend()
+        let M = 256
+        let fs = 256.0
+        let signal = Self.meanRemoved((0..<M).map { i in
+            let t = Double(i) / fs
+            return 2.0 * sin(2.0 * .pi * 8.0 * t) + 3.0 * sin(2.0 * .pi * 32.0 * t)
+        })
+
+        let psdSwift = pureSwift.powerSpectralDensity(signal, sampleRate: fs)
+        let psdAccel = accelerate.powerSpectralDensity(signal, sampleRate: fs)
+
+        #expect(psdSwift.count == psdAccel.count)
+        for (a, b) in zip(psdSwift, psdAccel) {
+            let tol = max(1e-9, max(abs(a), abs(b)) * 1e-9)
+            let delta = abs(a - b)
+            #expect(delta < tol)
+        }
+
+        // Both backends must satisfy Parseval, not just agree with each other.
+        let varSignal = Self.populationVariance(signal)
+        let intSwift = Self.parsevalIntegral(psd: psdSwift, sampleRate: fs)
+        let intAccel = Self.parsevalIntegral(psd: psdAccel, sampleRate: fs)
+        #expect(abs(intSwift - varSignal) / varSignal < 1e-12)
+        #expect(abs(intAccel - varSignal) / varSignal < 1e-12)
+    }
+    #endif
+
+    // MARK: - PSDBin convenience type
+
+    @Test("PSDBin: frequencies start at 0 and increment by Δf")
+    func psdBinFrequencySpacing() {
+        let backend = PureSwiftFFTBackend()
+        let signal = [Double](repeating: 0.0, count: 64)
+        let bins = backend.powerSpectralDensityBins(signal, sampleRate: 64.0)
+        let psd = backend.powerSpectralDensity(signal, sampleRate: 64.0)
+
+        #expect(bins.count == psd.count)
+        // Δf = fs / N_padded = 64 / 64 = 1.0
+        #expect(abs(bins[0].frequency - 0.0) < 1e-12)
+        #expect(abs(bins[1].frequency - 1.0) < 1e-12)
+        #expect(abs(bins[bins.count - 1].frequency - 32.0) < 1e-12)
+        // Power values match the bare PSD output
+        for (bin, p) in zip(bins, psd) {
+            #expect(bin.power == p)
+        }
+    }
+
+    @Test("PSDBin: empty input returns empty bins")
+    func psdBinEmptyInput() {
+        let backend = PureSwiftFFTBackend()
+        #expect(backend.powerSpectralDensityBins([], sampleRate: 100.0).isEmpty)
+    }
+}
diff --git a/Tests/Validation/PSD_Validation.swift b/Tests/Validation/PSD_Validation.swift
new file mode 100644
index 00000000..83f588f7
--- /dev/null
+++ b/Tests/Validation/PSD_Validation.swift
@@ -0,0 +1,258 @@
+// PSD_Validation.swift
+//
+// Standalone validation script for the PSD normalization upstream change
+// (v2.1.1). Hand-rolls the entire FFT + PSD pipeline with NO BusinessMath
+// dependency, then prints the values that the test suite will assert.
+//
+// Run with: swift Tests/Validation/PSD_Validation.swift
+//
+// The defining property of a correctly-normalized one-sided PSD is Parseval's
+// theorem:  σ²_time = Σ_k PSD[k] · Δf
+//
+// where Δf = fs / N_padded.
+//
+// This script verifies the formulas BEFORE any tests or production code touch
+// the package, so the test fixtures have an independent ground truth.
+
+import Foundation
+
+// MARK: - Naive radix-2 FFT (Cooley-Tukey, decimation in time)
+
+func nextPowerOf2(_ n: Int) -> 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..<bits {
+        r = (r << 1) | (v & 1)
+        v >>= 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..<n {
+        let j = bitReverse(i, bits: logN)
+        if i < j {
+            real.swapAt(i, j)
+            imag.swapAt(i, j)
+        }
+    }
+    var size = 2
+    while size <= n {
+        let half = size / 2
+        let angle = -2.0 * .pi / Double(size)
+        for start in stride(from: 0, to: n, by: size) {
+            for k in 0..<half {
+                let theta = angle * Double(k)
+                let c = cos(theta)
+                let s = sin(theta)
+                let evenIdx = start + k
+                let oddIdx = start + k + half
+                let tReal = c * real[oddIdx] - s * imag[oddIdx]
+                let tImag = s * real[oddIdx] + c * imag[oddIdx]
+                real[oddIdx] = real[evenIdx] - tReal
+                imag[oddIdx] = imag[evenIdx] - tImag
+                real[evenIdx] = real[evenIdx] + tReal
+                imag[evenIdx] = imag[evenIdx] + tImag
+            }
+        }
+        size *= 2
+    }
+}
+
+/// Raw power spectrum |X[k]|², matching BusinessMath's existing
+/// `powerSpectrum(_:)` semantics (zero-pads to next power of 2, returns
+/// N/2+1 bins).
+func rawPowerSpectrum(_ signal: [Double]) -> [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..<signal.count { real[i] = signal[i] }
+    fft(&real, &imag, n)
+    let bins = n / 2 + 1
+    var power = [Double](repeating: 0, count: bins)
+    for k in 0..<bins {
+        power[k] = real[k] * real[k] + imag[k] * imag[k]
+    }
+    return power
+}
+
+/// One-sided PSD with the M-vs-N normalization.
+/// PSD[k] = 2·|X[k]|² / (M·fs)   for typical bins
+/// PSD[0] = |X[0]|² / (M·fs)     (DC, not doubled)
+/// PSD[Nyq] = |X[Nyq]|² / (M·fs) (Nyquist, not doubled)
+func psd(_ signal: [Double], sampleRate: Double) -> [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..<nyquistBin {
+            out[k] = raw[k] * typicalFactor
+        }
+    }
+    return out
+}
+
+func variance(_ x: [Double]) -> 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..<M).map { i in A * sin(2.0 * .pi * f0 * Double(i) / fs) }
+    let timeVar = variance(x)
+    let p = psd(x, sampleRate: fs)
+    let integral = parsevalIntegral(psd: p, sampleRate: fs)
+    print("Time-domain variance:        \(timeVar)")
+    print("Analytic A²/2:                \(A * A / 2)")
+    print("PSD integral (should match):  \(integral)")
+    print("Relative error:               \(abs(integral - timeVar) / timeVar)")
+    print()
+}
+
+// MARK: Fixture 2 — Pure sine, M = 100 (NOT power of 2; padded to 128)
+do {
+    print("--- Fixture 2: Pure sine, M=100 (padded to 128), fs=100 Hz, f0=10 Hz ---")
+    let M = 100
+    let fs = 100.0
+    let f0 = 10.0
+    let A = 1.0
+    let x = (0..<M).map { i in A * sin(2.0 * .pi * f0 * Double(i) / fs) }
+    let timeVar = variance(x)
+    let p = psd(x, sampleRate: fs)
+    let integral = parsevalIntegral(psd: p, sampleRate: fs)
+    print("M (unpadded):                \(M)")
+    print("N (padded):                  \(nextPowerOf2(M))")
+    print("Time-domain variance:        \(timeVar)")
+    print("PSD integral:                \(integral)")
+    print("Relative error:              \(abs(integral - timeVar) / timeVar)")
+    print("Note: Some leakage expected because f0 doesn't land exactly on a bin")
+    print("after zero-padding (true bin freqs at fs/N_padded, not fs/M).")
+    print()
+}
+
+// MARK: Fixture 3 — Two-tone signal, M = 256
+do {
+    print("--- Fixture 3: Two tones, M=256, fs=256 Hz, f1=8 Hz A1=2, f2=32 Hz A2=3 ---")
+    let M = 256
+    let fs = 256.0
+    let f1 = 8.0
+    let f2 = 32.0
+    let A1 = 2.0
+    let A2 = 3.0
+    let x = (0..<M).map { i in
+        let t = Double(i) / fs
+        return A1 * sin(2.0 * .pi * f1 * t) + A2 * sin(2.0 * .pi * f2 * t)
+    }
+    let timeVar = variance(x)
+    let analyticVar = A1 * A1 / 2 + A2 * A2 / 2  // independent sines: variances add
+    let p = psd(x, sampleRate: fs)
+    let integral = parsevalIntegral(psd: p, sampleRate: fs)
+    print("Time-domain variance:        \(timeVar)")
+    print("Analytic A1²/2 + A2²/2:      \(analyticVar)")
+    print("PSD integral:                \(integral)")
+    print("Relative error:              \(abs(integral - timeVar) / timeVar)")
+    print()
+}
+
+// MARK: Fixture 4 — Zero-padding equivalence (with mean removal)
+do {
+    print("--- Fixture 4: Zero-padding equivalence (mean-removed signals) ---")
+    print("Same logical signal at M=64 (no padding) and M=50 (padded to 64)")
+    print("Mean is removed first, per the PSD method's contract.")
+    let fs = 64.0
+    let A = 1.5
+    let f0 = 8.0
+    func meanRemoved(_ x: [Double]) -> [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..<M).map { i in i.isMultiple(of: 2) ? 1.0 : -1.0 }
+    let timeVar = variance(x)  // 1.0
+    let p = psd(x, sampleRate: fs)
+    let integral = parsevalIntegral(psd: p, sampleRate: fs)
+    print("Time-domain variance:    \(timeVar)  (should be 1.0)")
+    print("PSD Nyquist bin:         \(p[p.count - 1])")
+    print("PSD other bins (max):    \(p.dropLast().max() ?? 0)")
+    print("PSD integral:            \(integral)  (should match variance)")
+    print("Relative error:          \(abs(integral - timeVar) / timeVar)")
+    print()
+}
+
+print("=== End validation ===")