diff --git a/Packages/PFFFT/Package.swift b/Packages/PFFFT/Package.swift index 6b05beb..051cf1a 100644 --- a/Packages/PFFFT/Package.swift +++ b/Packages/PFFFT/Package.swift @@ -29,7 +29,7 @@ let package = Package( ), .testTarget( name: "PFFFTTests", - dependencies: ["PFFFTLib"] + dependencies: ["PFFFT"] ), ] ) diff --git a/Packages/PFFFT/Sources/PFFFT/Buffer.swift b/Packages/PFFFT/Sources/PFFFT/Buffer.swift index d602ba2..2a2aaac 100644 --- a/Packages/PFFFT/Sources/PFFFT/Buffer.swift +++ b/Packages/PFFFT/Sources/PFFFT/Buffer.swift @@ -3,10 +3,12 @@ let bufferAlignment = 32 @frozen public struct Buffer: ~Copyable { let buffer: UnsafeMutableBufferPointer + var count: Int { buffer.count } + var baseAddress: UnsafeMutablePointer { buffer.baseAddress! } - init(unsafeUninitializedCapacity: Int) { + public init(capacity: Int) { buffer = UnsafeMutableRawBufferPointer.allocate( - byteCount: MemoryLayout.stride * unsafeUninitializedCapacity, + byteCount: MemoryLayout.stride * capacity, alignment: bufferAlignment ).bindMemory(to: Element.self) } @@ -15,15 +17,15 @@ public struct Buffer: ~Copyable { buffer.deallocate() } - func withUnsafeMutableBufferPointer(_ body: (UnsafeMutableBufferPointer) throws -> R) rethrows -> R { + public func withUnsafeMutableBufferPointer(_ body: (UnsafeMutableBufferPointer) throws -> R) rethrows -> R { try body(buffer) } - func withUnsafeBufferPointer(_ body: (UnsafeBufferPointer) throws -> R) rethrows -> R { + public func withUnsafeBufferPointer(_ body: (UnsafeBufferPointer) throws -> R) rethrows -> R { try body(UnsafeBufferPointer(buffer)) } - func withUnsafeMutableBytes(_ body: (UnsafeMutableRawBufferPointer) throws -> R) rethrows -> R { + public func withUnsafeMutableBytes(_ body: (UnsafeMutableRawBufferPointer) throws -> R) rethrows -> R { try body(UnsafeMutableRawBufferPointer(buffer)) } } diff --git a/Packages/PFFFT/Sources/PFFFT/FFTDoubleImpl.swift b/Packages/PFFFT/Sources/PFFFT/FFTDoubleImpl.swift index 172aefb..6d872f9 100644 --- a/Packages/PFFFT/Sources/PFFFT/FFTDoubleImpl.swift +++ b/Packages/PFFFT/Sources/PFFFT/FFTDoubleImpl.swift @@ -1,57 +1,69 @@ internal import PFFFTLib -struct FFTDoubleImpl: FFTImplProtocol { +public final class FFTDoubleImpl: FFTImplProtocol { let ptr: OpaquePointer + let n: Int + let type: FFTType - init?(n: Int, type: FFTType) { - guard let ptr = pffftd_new_setup(Int32(n), pffft_transform_t(type)) else { - return nil - } + public init(n: Int, type: FFTType) throws { + guard let ptr = pffftd_new_setup(Int32(n), pffft_transform_t(type)) else { throw FFTError.invalidSize } self.ptr = ptr + self.n = n + self.type = type } - func fft(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?, sign: FFTSign) { - let work: UnsafeMutablePointer! = switch work { - case let .some(b): b.buffer.baseAddress + public func fft(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?, sign: FFTSign) { + checkFftBufferCounts(n: n, type: type, input: input, output: output, work: work) + + let workAddress: UnsafeMutablePointer! = switch work { + case let .some(b): b.baseAddress case .none: nil } - - pffftd_transform_ordered(ptr, input.buffer.baseAddress, output.buffer.baseAddress, work, pffft_direction_t(sign)) + pffftd_transform_ordered(ptr, input.baseAddress, output.baseAddress, workAddress, pffft_direction_t(sign)) } - func fftUnordered(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?, sign: FFTSign) { - let work: UnsafeMutablePointer! = switch work { - case let .some(b): b.buffer.baseAddress + public func fftUnordered(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?, sign: FFTSign) { + checkFftBufferCounts(n: n, type: type, input: input, output: output, work: work) + + let workAddress: UnsafeMutablePointer! = switch work { + case let .some(b): b.baseAddress case .none: nil } - pffftd_transform(ptr, input.buffer.baseAddress, output.buffer.baseAddress, work, pffft_direction_t(sign)) + pffftd_transform(ptr, input.baseAddress, output.baseAddress, workAddress, pffft_direction_t(sign)) } - func zReorder(input: borrowing Buffer, output: borrowing Buffer, sign: FFTSign) { - pffftd_zreorder(ptr, input.buffer.baseAddress, output.buffer.baseAddress, pffft_direction_t(sign)) + public func zReorder(input: borrowing Buffer, output: borrowing Buffer, sign: FFTSign) { + checkFftBufferCounts(n: n, type: type, input: input, output: output, work: nil) + pffftd_zreorder(ptr, input.baseAddress, output.baseAddress, pffft_direction_t(sign)) } - func zConvolveAccumulate(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Double) { - pffftd_zconvolve_accumulate(ptr, dftA.buffer.baseAddress, dftB.buffer.baseAddress, dftAB.buffer.baseAddress, scaling) + public func zConvolveAccumulate(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Double) { + checkConvolveBufferCounts(n: n, type: type, a: dftA, b: dftB, ab: dftAB) + pffftd_zconvolve_accumulate(ptr, dftA.baseAddress, dftB.baseAddress, dftAB.baseAddress, scaling) } - func zConvolve(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Double) { - pffftd_zconvolve_no_accu(ptr, dftA.buffer.baseAddress, dftB.buffer.baseAddress, dftAB.buffer.baseAddress, scaling) + public func zConvolve(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Double) { + checkConvolveBufferCounts(n: n, type: type, a: dftA, b: dftB, ab: dftAB) + pffftd_zconvolve_no_accu(ptr, dftA.baseAddress, dftB.baseAddress, dftAB.baseAddress, scaling) } - func minFftSize(for type: FFTType) -> Int { + public static func minFftSize(for type: FFTType) -> Int { Int(pffftd_min_fft_size(pffft_transform_t(type))) } - static func isValidSize(_ n: Int, for type: FFTType) -> Bool { + public static func isValidSize(_ n: Int, for type: FFTType) -> Bool { pffftd_is_valid_size(Int32(n), pffft_transform_t(type)) != 0 } - static func nearestValidSize(_ n: Int, for type: FFTType, higher: Bool) -> Int { + public static func nearestValidSize(_ n: Int, for type: FFTType, higher: Bool) -> Int { Int(pffftd_nearest_transform_size(Int32(n), pffft_transform_t(type), higher ? 1 : 0)) } - static func simdArch() -> String { - String(cString: pffftd_simd_arch()) + deinit { + pffftd_destroy_setup(ptr) } } + +extension Double: FFTElement { + public typealias FFTImpl = FFTDoubleImpl +} diff --git a/Packages/PFFFT/Sources/PFFFT/FFTFloatImpl.swift b/Packages/PFFFT/Sources/PFFFT/FFTFloatImpl.swift index ef6d897..c61008f 100644 --- a/Packages/PFFFT/Sources/PFFFT/FFTFloatImpl.swift +++ b/Packages/PFFFT/Sources/PFFFT/FFTFloatImpl.swift @@ -1,57 +1,69 @@ internal import PFFFTLib -struct FFTFloatImpl: FFTImplProtocol { +public final class FFTFloatImpl: FFTImplProtocol { let ptr: OpaquePointer + let n: Int + let type: FFTType - init?(n: Int, type: FFTType) { - guard let ptr = pffft_new_setup(Int32(n), pffft_transform_t(type)) else { - return nil - } + public init(n: Int, type: FFTType) throws { + guard let ptr = pffft_new_setup(Int32(n), pffft_transform_t(type)) else { throw FFTError.invalidSize } self.ptr = ptr + self.n = n + self.type = type } - func fft(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?, sign: FFTSign) { - let work: UnsafeMutablePointer! = switch work { - case let .some(b): b.buffer.baseAddress + public func fft(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?, sign: FFTSign) { + checkFftBufferCounts(n: n, type: type, input: input, output: output, work: work) + + let workAddress: UnsafeMutablePointer! = switch work { + case let .some(b): b.baseAddress case .none: nil } - - pffft_transform_ordered(ptr, input.buffer.baseAddress, output.buffer.baseAddress, work, pffft_direction_t(sign)) + pffft_transform_ordered(ptr, input.baseAddress, output.baseAddress, workAddress, pffft_direction_t(sign)) } - func fftUnordered(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?, sign: FFTSign) { - let work: UnsafeMutablePointer! = switch work { - case let .some(b): b.buffer.baseAddress + public func fftUnordered(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?, sign: FFTSign) { + checkFftBufferCounts(n: n, type: type, input: input, output: output, work: work) + + let workAddress: UnsafeMutablePointer! = switch work { + case let .some(b): b.baseAddress case .none: nil } - pffft_transform(ptr, input.buffer.baseAddress, output.buffer.baseAddress, work, pffft_direction_t(sign)) + pffft_transform(ptr, input.baseAddress, output.baseAddress, workAddress, pffft_direction_t(sign)) } - func zReorder(input: borrowing Buffer, output: borrowing Buffer, sign: FFTSign) { - pffft_zreorder(ptr, input.buffer.baseAddress, output.buffer.baseAddress, pffft_direction_t(sign)) + public func zReorder(input: borrowing Buffer, output: borrowing Buffer, sign: FFTSign) { + checkFftBufferCounts(n: n, type: type, input: input, output: output, work: nil) + pffft_zreorder(ptr, input.baseAddress, output.baseAddress, pffft_direction_t(sign)) } - func zConvolveAccumulate(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Float) { - pffft_zconvolve_accumulate(ptr, dftA.buffer.baseAddress, dftB.buffer.baseAddress, dftAB.buffer.baseAddress, scaling) + public func zConvolveAccumulate(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Float) { + checkConvolveBufferCounts(n: n, type: type, a: dftA, b: dftB, ab: dftAB) + pffft_zconvolve_accumulate(ptr, dftA.baseAddress, dftB.baseAddress, dftAB.baseAddress, scaling) } - func zConvolve(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Float) { - pffft_zconvolve_no_accu(ptr, dftA.buffer.baseAddress, dftB.buffer.baseAddress, dftAB.buffer.baseAddress, scaling) + public func zConvolve(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Float) { + checkConvolveBufferCounts(n: n, type: type, a: dftA, b: dftB, ab: dftAB) + pffft_zconvolve_no_accu(ptr, dftA.baseAddress, dftB.baseAddress, dftAB.baseAddress, scaling) } - func minFftSize(for type: FFTType) -> Int { + public static func minFftSize(for type: FFTType) -> Int { Int(pffft_min_fft_size(pffft_transform_t(type))) } - static func isValidSize(_ n: Int, for type: FFTType) -> Bool { + public static func isValidSize(_ n: Int, for type: FFTType) -> Bool { pffft_is_valid_size(Int32(n), pffft_transform_t(type)) != 0 } - static func nearestValidSize(_ n: Int, for type: FFTType, higher: Bool) -> Int { + public static func nearestValidSize(_ n: Int, for type: FFTType, higher: Bool) -> Int { Int(pffft_nearest_transform_size(Int32(n), pffft_transform_t(type), higher ? 1 : 0)) } - static func simdArch() -> String { - String(cString: pffft_simd_arch()) + deinit { + pffft_destroy_setup(ptr) } } + +extension Float: FFTElement { + public typealias FFTImpl = FFTFloatImpl +} diff --git a/Packages/PFFFT/Sources/PFFFT/PFFFT.swift b/Packages/PFFFT/Sources/PFFFT/PFFFT.swift index 87688ab..002a7f6 100644 --- a/Packages/PFFFT/Sources/PFFFT/PFFFT.swift +++ b/Packages/PFFFT/Sources/PFFFT/PFFFT.swift @@ -1,29 +1,111 @@ internal import PFFFTLib -enum FFTType { +@frozen +public enum FFTType { case real case complex } -enum FFTSign: Int { +@frozen +public enum FFTSign: Int { case forward = -1 case backward = 1 } -protocol FFTImplProtocol { +public enum FFTError: Error { + case invalidSize +} + +public protocol FFTElement { + associatedtype FFTImpl: FFTImplProtocol +} + +public protocol FFTImplProtocol: ~Copyable { associatedtype Element - init?(n: Int, type: FFTType) + /// Initialize the FFT implementation with the given size and type. + /// - Parameters: + /// - n: The size of the FFT. + /// - type: The type of FFT. + /// - Throws: `FFTError.invalidSize` if the size is invalid. + init(n: Int, type: FFTType) throws + /// Perform a forward or backward FFT on the input buffer. + /// + /// The input and output buffers may be the same. + /// The data is stores in order as expected (interleaved complex components ordered by frequency). + /// The input and output buffer must have a capacity of at least `n` for real FFTs and `2 * n` for complex FFTs. + /// A fatal error will occur if any buffer is too small. + /// + /// For a real forward transform with real input, the output array is organized as follows: + /// index k > 2 where k is even is the real part of the k/2-th complex coefficient. + /// index k > 2 where k is odd is the imaginary part of the k/2-th complex coefficient. + /// index k = 0 is the real part of the 0 frequency (DC) coefficient. + /// index k = 1 is the real part of the Nyquist coefficient. + /// + /// Transforms are not scaled. fft_backward(fft_forward(x)) == n * x. + /// + /// - Parameters: + /// - input: The input buffer. + /// - output: The output buffer. + /// - work: An optional work buffer. Must have capacity of at least `n` for real FFTs and `2 * n` for complex FFTs. + /// - sign: The direction of the FFT. func fft(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?, sign: FFTSign) + + /// Perform a forward or backward FFT on the input buffer, with implementation defined order. + /// + /// This function behaves similarly to `fft` however the z-domain data is stored in most efficient ordering, + /// which is suitable for transforming back with this function, or for convolution. + /// - Parameters: + /// - input: The input buffer. + /// - output: The output buffer. + /// - work: An optional work buffer. Must have capacity of at least `n` for real FFTs and `2 * n` for complex FFTs. + /// - sign: The direction of the FFT. func fftUnordered(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?, sign: FFTSign) + func zReorder(input: borrowing Buffer, output: borrowing Buffer, sign: FFTSign) + + /// Perform a convolution of two complex signals in the frequency domain. + /// + /// Multiplies frequency domain components of `dftA` and `dftB` and accumulates the result in `dftAB`. + /// The operation performed is `dftAB += (dftA * dftB) * scaling`. + /// - Parameters: + /// - dftA: The first input buffer of frequency domain data. + /// - dftB: The second input buffer of frequency domain data. + /// - dftAB: The output buffer of frequency domain data. + /// - scaling: The scaling factor to apply to the result. func zConvolveAccumulate(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Element) + + /// Perform a convolution of two complex signals in the frequency domain. + /// + /// Multiplies frequency domain components of `dftA` and `dftB` and stores the result in `dftAB`. + /// The operation performed is `dftAB = (dftA * dftB) * scaling`. + /// - Parameters: + /// - dftA: The first input buffer of frequency domain data. + /// - dftB: The second input buffer of frequency domain data. + /// - dftAB: The output buffer of frequency domain data. + /// - scaling: The scaling factor to apply to the result. func zConvolve(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Element) - func minFftSize(for type: FFTType) -> Int + + /// Returns the minimum FFT size for the given type. + /// + /// - Parameter type: The type of FFT. + static func minFftSize(for type: FFTType) -> Int + + /// Returns whether the given size is valid for the given type. + /// + /// The PFFFT library requires `n` to be factorizable to `minFftSize` with factors of 2, 3, 5. + /// - Parameters: + /// - n: The size to check. + /// - type: The type of FFT. static func isValidSize(_ n: Int, for type: FFTType) -> Bool + + /// Returns the nearest valid size for the given type. static func nearestValidSize(_ n: Int, for type: FFTType, higher: Bool) -> Int - static func simdArch() -> String +} + +public var simdArch: String { + String(cString: pffft_simd_arch()) } extension pffft_transform_t { @@ -45,3 +127,38 @@ extension pffft_direction_t { } } } + +@inline(__always) +func checkFftBufferCounts(n: Int, type: FFTType, input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?) { + let minSize = type == .real ? n : 2 * n + + let work: UnsafeMutableBufferPointer! = switch work { + case let .some(b): b.buffer + case .none: nil + } + + guard input.count >= minSize else { + fatalError("input buffer too small") + } + guard output.count >= minSize else { + fatalError("output buffer too small") + } + guard work == nil || work.count >= minSize else { + fatalError("work buffer too small") + } +} + +@inline(__always) +func checkConvolveBufferCounts(n: Int, type: FFTType, a: borrowing Buffer, b: borrowing Buffer, ab: borrowing Buffer) { + let minSize = type == .real ? n : 2 * n + + guard a.count >= minSize else { + fatalError("a buffer too small") + } + guard b.count >= minSize else { + fatalError("b buffer too small") + } + guard ab.count >= minSize else { + fatalError("ab buffer too small") + } +} diff --git a/Packages/PFFFT/Sources/PFFFT/SetupCache.swift b/Packages/PFFFT/Sources/PFFFT/SetupCache.swift index e45afb3..10bbc29 100644 --- a/Packages/PFFFT/Sources/PFFFT/SetupCache.swift +++ b/Packages/PFFFT/Sources/PFFFT/SetupCache.swift @@ -5,13 +5,15 @@ struct CacheKey : Hashable { let type: FFTType } -class SetupCache : @unchecked Sendable { - private var cache: [CacheKey: Impl?] = [:] +public class SetupCache : @unchecked Sendable { + private var cache: [CacheKey: Element.FFTImpl?] = [:] private let queue = DispatchQueue(label: String(describing: SetupCache.self), attributes: .concurrent) - func get(for n: Int, type: FFTType) -> Impl? { - var setup: Impl?? + public init() {} + + public func get(for n: Int, type: FFTType) throws -> Element.FFTImpl { + var setup: Element.FFTImpl?? queue.sync { setup = cache[CacheKey(n: n, type: type)] } @@ -19,12 +21,19 @@ class SetupCache : @unchecked Sendable { queue.sync(flags: .barrier) { setup = cache[CacheKey(n: n, type: type)] if setup == nil { - setup = Impl(n: n, type: type) - cache[CacheKey(n: n, type: type)] = setup + let entry = try? Element.FFTImpl(n: n, type: type) + cache[CacheKey(n: n, type: type)] = entry + setup = entry } } } - return setup! + guard let s = setup! else { throw FFTError.invalidSize } + return s } + public func clear() { + queue.sync(flags: .barrier) { + cache.removeAll() + } + } } diff --git a/Packages/PFFFT/Tests/PFFFTTests/PFFFTTests.swift b/Packages/PFFFT/Tests/PFFFTTests/PFFFTTests.swift index 1ca34d3..5a4021e 100644 --- a/Packages/PFFFT/Tests/PFFFTTests/PFFFTTests.swift +++ b/Packages/PFFFT/Tests/PFFFTTests/PFFFTTests.swift @@ -2,5 +2,14 @@ import Testing @testable import PFFFT @Test func example() async throws { + let cache = PFFFT.SetupCache() + + let ff = try cache.get(for: 1024, type: .real) + + let input = Buffer(capacity: 1023) + let output = Buffer(capacity: 1024) + + ff.fft(input: input, output: output, work: nil, sign: .forward) + // Write your test here and use APIs like `#expect(...)` to check expected conditions. } diff --git a/Sources/EcgSynKit/EcgSynKit.swift b/Sources/EcgSynKit/EcgSynKit.swift index a3c1b29..81e2ab7 100644 --- a/Sources/EcgSynKit/EcgSynKit.swift +++ b/Sources/EcgSynKit/EcgSynKit.swift @@ -1,6 +1,7 @@ import Algorithms import Foundation import PFFFTLib +import PFFFT struct Parameters { /// The number of beats to simulate. @@ -52,39 +53,6 @@ struct Parameters { let fhistd = 0.01 } -class PfftSetupCache: @unchecked Sendable { - private var cache: [Int: OpaquePointer?] = [:] - - private let queue = DispatchQueue(label: String(describing: PfftSetupCache.self), attributes: .concurrent) - - func get(for nrr: Int) -> OpaquePointer? { - var setup: OpaquePointer?? - queue.sync { - setup = cache[nrr] - } - if setup == nil { - queue.sync(flags: .barrier) { - setup = cache[nrr] - if setup == nil { - setup = pffftd_new_setup(Int32(nrr), PFFFT_REAL) - cache[nrr] = setup - } - } - } - return setup! - } - - deinit { - for (_, setup) in cache { - if setup != nil { - pffftd_destroy_setup(setup) - } - } - } - - static let shared = PfftSetupCache() -} - struct EcgDerive { let rrpc: [Double] } @@ -95,6 +63,8 @@ func stdev(_ data: [Double]) -> Double { return sqrt(data.lazy.map { ($0 - mean) * ($0 - mean) }.reduce(0.0, +) / (n - 1)) } +let fftlib = SetupCache() + func rrprocess(params: Parameters, nrr: Int) -> [Double] { let w1 = 2.0 * .pi * params.flo let w2 = 2.0 * .pi * params.fhi @@ -110,11 +80,11 @@ func rrprocess(params: Parameters, nrr: Int) -> [Double] { let sf = Double(params.sfInternal) let df = sf / Double(nrr) - let fft = PfftSetupCache.shared.get(for: nrr) - var rr = withUnsafeTemporaryAllocation(byteCount: (nrr + 2) * MemoryLayout.stride, alignment: 32) { - let swc = $0.bindMemory(to: Double.self) + let input = Buffer(capacity: (nrr + 2)) + let output = Buffer(capacity: nrr) + input.withUnsafeMutableBufferPointer { swc in for i in 0 ..< nrr / 2 + 1 { let w = df * Double(i) * 2.0 * .pi let dw1 = w - w1 @@ -131,13 +101,12 @@ func rrprocess(params: Parameters, nrr: Int) -> [Double] { // pack Nyquist frequency real to imaginary of DC swc[1] = swc[nrr] + } - return withUnsafeTemporaryAllocation(byteCount: nrr * MemoryLayout.stride, alignment: 32) { - let outptr = $0.bindMemory(to: Double.self) + try! fftlib.get(for: nrr, type: .real).fft(input: input, output: output, work: nil, sign: .backward) - pffftd_transform_ordered(fft, swc.baseAddress, outptr.baseAddress, nil, PFFFT_BACKWARD) - return outptr.map { $0 * 1.0 / Double(nrr) } - } + var rr = output.withUnsafeMutableBufferPointer { outptr in + return outptr.map { $0 * 1.0 / Double(nrr) } } let xstd = stdev(rr)