diff --git a/Package.swift b/Package.swift index 030afe0..f1f6896 100644 --- a/Package.swift +++ b/Package.swift @@ -14,7 +14,6 @@ let package = Package( ], dependencies: [ .package(url: "https://github.com/apple/swift-algorithms", from: "1.2.0"), - .package(path: "./Packages/KissFFT/"), .package(path: "./Packages/PFFFT/"), ], targets: [ diff --git a/Packages/PFFFT/Package.swift b/Packages/PFFFT/Package.swift index 051cf1a..dce638b 100644 --- a/Packages/PFFFT/Package.swift +++ b/Packages/PFFFT/Package.swift @@ -6,12 +6,14 @@ import PackageDescription let package = Package( name: "PFFFT", products: [ - // Products define the executables and libraries a package produces, making them visible to other packages. .library( name: "PFFFT", targets: ["PFFFT", "PFFFTLib"] ), ], + dependencies: [ + .package(url: "https://github.com/apple/swift-numerics", from: "1.0.0"), + ], targets: [ .target( name: "PFFFTLib", @@ -25,7 +27,7 @@ let package = Package( ), .target( name: "PFFFT", - dependencies: ["PFFFTLib"] + dependencies: ["PFFFTLib", .product(name: "Numerics", package: "swift-numerics")] ), .testTarget( name: "PFFFTTests", diff --git a/Packages/PFFFT/Sources/PFFFT/Buffer.swift b/Packages/PFFFT/Sources/PFFFT/Buffer.swift index 2a2aaac..67c1988 100644 --- a/Packages/PFFFT/Sources/PFFFT/Buffer.swift +++ b/Packages/PFFFT/Sources/PFFFT/Buffer.swift @@ -1,27 +1,27 @@ let bufferAlignment = 32 @frozen -public struct Buffer: ~Copyable { - let buffer: UnsafeMutableBufferPointer +public struct Buffer: ~Copyable { + let buffer: UnsafeMutableBufferPointer var count: Int { buffer.count } - var baseAddress: UnsafeMutablePointer { buffer.baseAddress! } + var baseAddress: UnsafeMutablePointer { buffer.baseAddress! } public init(capacity: Int) { buffer = UnsafeMutableRawBufferPointer.allocate( - byteCount: MemoryLayout.stride * capacity, + byteCount: MemoryLayout.stride * capacity, alignment: bufferAlignment - ).bindMemory(to: Element.self) + ).bindMemory(to: T.self) } deinit { buffer.deallocate() } - public func withUnsafeMutableBufferPointer(_ body: (UnsafeMutableBufferPointer) throws -> R) rethrows -> R { + public func withUnsafeMutableBufferPointer(_ body: (UnsafeMutableBufferPointer) throws -> R) rethrows -> R { try body(buffer) } - public func withUnsafeBufferPointer(_ body: (UnsafeBufferPointer) throws -> R) rethrows -> R { + public func withUnsafeBufferPointer(_ body: (UnsafeBufferPointer) throws -> R) rethrows -> R { try body(UnsafeBufferPointer(buffer)) } diff --git a/Packages/PFFFT/Sources/PFFFT/FFT.swift b/Packages/PFFFT/Sources/PFFFT/FFT.swift new file mode 100644 index 0000000..c832bea --- /dev/null +++ b/Packages/PFFFT/Sources/PFFFT/FFT.swift @@ -0,0 +1,276 @@ +internal import PFFFTLib +import ComplexModule +import RealModule + +@frozen +public enum FFTType { + case real + case complex +} + +@frozen +public enum FFTSign: Int { + case forward = -1 + case backward = 1 +} + +public enum FFTError: Error { + case invalidSize +} + +public protocol FFTElemental { + associatedtype ScalarType: FFTScalar + associatedtype ComplexType = Complex + + static func setupPfft(_ n: Int, _ type: FFTType) throws -> OpaquePointer + static func pffftMinFftSize(_ type: FFTType) -> Int + static func pffftIsValidSize(_ n: Int, _ type: FFTType) -> Bool + static func pffftNearestValidSize(_ n: Int, _ type: FFTType, _ higher: Bool) -> Int +} + +public protocol FFTScalar: Real { + static func pffftTransformOrdered(_ ptr: OpaquePointer, _ input: UnsafeMutablePointer, _ output: UnsafeMutablePointer, _ work: UnsafeMutablePointer?, _ dir: FFTSign) + static func pffftTransform(_ ptr: OpaquePointer, _ input: UnsafeMutablePointer, _ output: UnsafeMutablePointer, _ work: UnsafeMutablePointer?, _ dir: FFTSign) + static func pffftZreorder(_ ptr: OpaquePointer, _ input: UnsafeMutablePointer, _ output: UnsafeMutablePointer, _ dir: FFTSign) + static func pffftZconvolveAccumulate(_ ptr: OpaquePointer, _ dftA: UnsafeMutablePointer, _ dftB: UnsafeMutablePointer, _ dftAB: UnsafeMutablePointer, _ scaling: Self) + static func pffftZconvolveNoAccu(_ ptr: OpaquePointer, _ dftA: UnsafeMutablePointer, _ dftB: UnsafeMutablePointer, _ dftAB: UnsafeMutablePointer, _ scaling: Self) +} + +extension Float: FFTScalar { + public static func pffftTransformOrdered(_ ptr: OpaquePointer, _ input: UnsafeMutablePointer, _ output: UnsafeMutablePointer, _ work: UnsafeMutablePointer?, _ dir: FFTSign) { + pffft_transform_ordered(ptr, input, output, work, pffft_direction_t(dir)) + } + + public static func pffftTransform(_ ptr: OpaquePointer, _ input: UnsafeMutablePointer, _ output: UnsafeMutablePointer, _ work: UnsafeMutablePointer?, _ dir: FFTSign) { + pffft_transform(ptr, input, output, work, pffft_direction_t(dir)) + } + + public static func pffftZreorder(_ ptr: OpaquePointer, _ input: UnsafeMutablePointer, _ output: UnsafeMutablePointer, _ dir: FFTSign) { + pffft_zreorder(ptr, input, output, pffft_direction_t(dir)) + } + + public static func pffftZconvolveAccumulate(_ ptr: OpaquePointer, _ dftA: UnsafeMutablePointer, _ dftB: UnsafeMutablePointer, _ dftAB: UnsafeMutablePointer, _ scaling: Self) { + pffft_zconvolve_accumulate(ptr, dftA, dftB, dftAB, scaling) + } + + public static func pffftZconvolveNoAccu(_ ptr: OpaquePointer, _ dftA: UnsafeMutablePointer, _ dftB: UnsafeMutablePointer, _ dftAB: UnsafeMutablePointer, _ scaling: Self) { + pffft_zconvolve_no_accu(ptr, dftA, dftB, dftAB, scaling) + } +} + +extension Double: FFTScalar { + public static func pffftTransformOrdered(_ ptr: OpaquePointer, _ input: UnsafeMutablePointer, _ output: UnsafeMutablePointer, _ work: UnsafeMutablePointer?, _ dir: FFTSign) { + pffftd_transform_ordered(ptr, input, output, work, pffft_direction_t(dir)) + } + + public static func pffftTransform(_ ptr: OpaquePointer, _ input: UnsafeMutablePointer, _ output: UnsafeMutablePointer, _ work: UnsafeMutablePointer?, _ dir: FFTSign) { + pffftd_transform(ptr, input, output, work, pffft_direction_t(dir)) + } + + public static func pffftZreorder(_ ptr: OpaquePointer, _ input: UnsafeMutablePointer, _ output: UnsafeMutablePointer, _ dir: FFTSign) { + pffftd_zreorder(ptr, input, output, pffft_direction_t(dir)) + } + + public static func pffftZconvolveAccumulate(_ ptr: OpaquePointer, _ dftA: UnsafeMutablePointer, _ dftB: UnsafeMutablePointer, _ dftAB: UnsafeMutablePointer, _ scaling: Double) { + pffftd_zconvolve_accumulate(ptr, dftA, dftB, dftAB, scaling) + } + + public static func pffftZconvolveNoAccu(_ ptr: OpaquePointer, _ dftA: UnsafeMutablePointer, _ dftB: UnsafeMutablePointer, _ dftAB: UnsafeMutablePointer, _ scaling: Double) { + pffftd_zconvolve_no_accu(ptr, dftA, dftB, dftAB, scaling) + } +} + +extension Complex: FFTElemental where RealType: FFTElemental & FFTScalar { + public typealias ScalarType = RealType + + public static func setupPfft(_ n: Int, _: FFTType) throws -> OpaquePointer { + return try ScalarType.setupPfft(n, .complex) + } + + public static func pffftMinFftSize(_: FFTType) -> Int { + return ScalarType.pffftMinFftSize(.complex) + } + + public static func pffftIsValidSize(_ n: Int, _: FFTType) -> Bool { + return ScalarType.pffftIsValidSize(n, .complex) + } + + public static func pffftNearestValidSize(_ n: Int, _: FFTType, _ higher: Bool) -> Int { + return ScalarType.pffftNearestValidSize(n, .complex, higher) + } +} + +extension Double: FFTElemental { + public typealias ScalarType = Double + + public static func setupPfft(_ n: Int, _ type: FFTType) throws -> OpaquePointer { + guard let ptr = pffftd_new_setup(Int32(n), pffft_transform_t(type)) else { throw FFTError.invalidSize } + return ptr + } + + public static func pffftMinFftSize(_ type: FFTType) -> Int { + Int(pffftd_min_fft_size(pffft_transform_t(type))) + } + + public static func pffftIsValidSize(_ n: Int, _ type: FFTType) -> Bool { + pffftd_is_valid_size(Int32(n), pffft_transform_t(type)) != 0 + } + + public static func pffftNearestValidSize(_ n: Int, _ type: FFTType, _ higher: Bool) -> Int { + Int(pffftd_nearest_transform_size(Int32(n), pffft_transform_t(type), higher ? 1 : 0)) + } +} + +extension Float: FFTElemental { + public typealias ScalarType = Float + + public static func setupPfft(_ n: Int, _ type: FFTType) throws -> OpaquePointer { + guard let ptr = pffft_new_setup(Int32(n), pffft_transform_t(type)) else { throw FFTError.invalidSize } + return ptr + } + + public static func pffftMinFftSize(_ type: FFTType) -> Int { + Int(pffft_min_fft_size(pffft_transform_t(type))) + } + + public static func pffftIsValidSize(_ n: Int, _ type: FFTType) -> Bool { + pffft_is_valid_size(Int32(n), pffft_transform_t(type)) != 0 + } + + public static func pffftNearestValidSize(_ n: Int, _ type: FFTType, _ higher: Bool) -> Int { + Int(pffft_nearest_transform_size(Int32(n), pffft_transform_t(type), higher ? 1 : 0)) + } +} + +public final class FFT { + public typealias ComplexType = T.ComplexType + public typealias ScalarType = T.ScalarType + + let ptr: OpaquePointer + let n: Int + let work: Buffer? + + /// 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. + public init(n: Int) throws { + ptr = try T.setupPfft(n, .real) + self.n = n + let capacity = T.self == Complex.self ? 2 * n : n + work = n > 4096 ? Buffer(capacity: capacity) : nil + } + + public func makeSignalBuffer(extra _: Int = 0) -> Buffer { + Buffer(capacity: n) + } + + public func makeSpectrumBuffer(extra _: Int = 0) -> Buffer { + Buffer(capacity: n) + } + + @inline(__always) + func toAddress(_ work: borrowing Buffer?) -> UnsafeMutablePointer? { + switch work { + case let .some(b): return b.baseAddress + case .none: return nil + } + } + + @inline(__always) + func rebind(_ buffer: borrowing Buffer) -> UnsafeMutablePointer! { + UnsafeMutableRawBufferPointer(buffer.buffer).bindMemory(to: ScalarType.self).baseAddress + } + + @inline(__always) + func checkFftBufferCounts(signal: borrowing Buffer, spectrum: borrowing Buffer) { + guard signal.count >= n else { + fatalError("signal buffer too small") + } + guard spectrum.count >= n else { + fatalError("spectrum buffer too small") + } + } + + @inline(__always) + func checkFftInternalLayoutBufferCounts(signal: borrowing Buffer, spectrum: borrowing Buffer) { + guard signal.count >= n else { + fatalError("signal buffer too small") + } + guard spectrum.count >= (T.self == Complex.self ? 2 * n : n) else { + fatalError("spectrum buffer too small") + } + } + + @inline(__always) + func checkConvolveBufferCounts(a: borrowing Buffer, b: borrowing Buffer, ab: borrowing Buffer) { + let minCount = T.self == Complex.self ? 2 * n : n + + guard a.count >= minCount else { + fatalError("a buffer too small") + } + guard b.count >= minCount else { + fatalError("b buffer too small") + } + guard ab.count >= minCount else { + fatalError("ab buffer too small") + } + } + + public func forward(signal: borrowing Buffer, spectrum: borrowing Buffer) { + checkFftBufferCounts(signal: signal, spectrum: spectrum) + ScalarType.pffftTransformOrdered(ptr, rebind(signal), rebind(spectrum), toAddress(work), .forward) + } + + public func inverse(spectrum: borrowing Buffer, signal: borrowing Buffer) { + checkFftBufferCounts(signal: signal, spectrum: spectrum) + ScalarType.pffftTransformOrdered(ptr, rebind(spectrum), rebind(signal), toAddress(work), .backward) + } + + public func forwardToInternalLayout(signal: borrowing Buffer, spectrum: borrowing Buffer) { + checkFftInternalLayoutBufferCounts(signal: signal, spectrum: spectrum) + ScalarType.pffftTransform(ptr, rebind(signal), spectrum.baseAddress, toAddress(work), .forward) + } + + public func inverseFromInternalLayout(spectrum: borrowing Buffer, signal: borrowing Buffer) { + checkFftInternalLayoutBufferCounts(signal: signal, spectrum: spectrum) + ScalarType.pffftTransform(ptr, spectrum.baseAddress, rebind(signal), toAddress(work), .backward) + } + + public func reorder(spectrum: borrowing Buffer, output: borrowing Buffer) { + guard spectrum.count >= (T.self == Complex.self ? 2 * n : n) else { + fatalError("signal buffer too small") + } + guard output.count >= n else { + fatalError("outputbuffer too small") + } + ScalarType.pffftZreorder(ptr, spectrum.baseAddress, rebind(output), .forward) + } + + public func convolve(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: ScalarType) { + checkConvolveBufferCounts(a: dftA, b: dftB, ab: dftAB) + ScalarType.pffftZconvolveNoAccu(ptr, dftA.baseAddress, dftB.baseAddress, dftAB.baseAddress, scaling) + } + + public func convolveAccumulate(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: ScalarType) { + checkConvolveBufferCounts(a: dftA, b: dftB, ab: dftAB) + ScalarType.pffftZconvolveAccumulate(ptr, dftA.baseAddress, dftB.baseAddress, dftAB.baseAddress, scaling) + } + + public static func minFftSize() -> Int { + T.pffftMinFftSize(.real) + } + + public static func isValidSize(_ n: Int) -> Bool { + T.pffftIsValidSize(n, .real) + } + + public static func nearestValidSize(_ n: Int, higher: Bool) -> Int { + T.pffftNearestValidSize(n, .real, higher) + } + + deinit { + pffftd_destroy_setup(ptr) + } +} diff --git a/Packages/PFFFT/Sources/PFFFT/FFTDouble.swift b/Packages/PFFFT/Sources/PFFFT/FFTDouble.swift deleted file mode 100644 index 7b2699e..0000000 --- a/Packages/PFFFT/Sources/PFFFT/FFTDouble.swift +++ /dev/null @@ -1,94 +0,0 @@ -internal import PFFFTLib - -public final class FFTDouble: PFFFTProtocol { - let ptr: OpaquePointer - let n: Int - let type: FFTType - - static let sharedCache = SetupCache() - - public static func setup(for n: Int, type: FFTType) throws -> FFTDouble { - try sharedCache.get(for: n, type: type) - } - - 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 transform(_ input: borrowing Buffer, _ output: borrowing Buffer, _ work: borrowing Buffer?, _ dir: pffft_direction_t) { - 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.baseAddress, output.baseAddress, workAddress, dir) - } - - func transformUnordered(_ input: borrowing Buffer, _ output: borrowing Buffer, _ work: borrowing Buffer?, _ dir: pffft_direction_t) { - 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.baseAddress, output.baseAddress, workAddress, dir) - } - - - public func forward(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?) { - transform(input, output, work, PFFFT_FORWARD) - } - - public func inverse(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?) { - transform(input, output, work, PFFFT_BACKWARD) - } - - public func forwardUnordered(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?) { - transformUnordered(input, output, work, PFFFT_FORWARD) - } - - public func inverseUnordered(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?) { - transformUnordered(input, output, work, PFFFT_BACKWARD) - } - - - public func reorderSpectrum(input: borrowing Buffer, output: borrowing Buffer) { - checkFftBufferCounts(n: n, type: type, input: input, output: output, work: nil) - pffftd_zreorder(ptr, input.baseAddress, output.baseAddress, PFFFT_FORWARD) - } - - public func convolveAccumulate(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) - } - - public func convolve(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) - } - - public static func minFftSize(for type: FFTType) -> Int { - Int(pffftd_min_fft_size(pffft_transform_t(type))) - } - - public static func isValidSize(_ n: Int, for type: FFTType) -> Bool { - pffftd_is_valid_size(Int32(n), pffft_transform_t(type)) != 0 - } - - 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)) - } - - deinit { - pffftd_destroy_setup(ptr) - } -} - -extension Double: FFTElement { - public typealias FFTImpl = FFTDouble -} diff --git a/Packages/PFFFT/Sources/PFFFT/FFTFloat.swift b/Packages/PFFFT/Sources/PFFFT/FFTFloat.swift deleted file mode 100644 index 9b0d05a..0000000 --- a/Packages/PFFFT/Sources/PFFFT/FFTFloat.swift +++ /dev/null @@ -1,75 +0,0 @@ -internal import PFFFTLib - -public final class FFTFloat: PFFFTProtocol { - let ptr: OpaquePointer - let n: Int - let type: FFTType - - static let sharedCache = SetupCache() - - public static func setup(for n: Int, type: FFTType) throws -> FFTFloat { - try sharedCache.get(for: n, type: type) - } - - 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 - } - - public func forward(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.baseAddress, output.baseAddress, workAddress, pffft_direction_t(sign)) - } - - 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.baseAddress, output.baseAddress, workAddress, pffft_direction_t(sign)) - } - - public func reorderSpectrum(input: borrowing Buffer, output: borrowing Buffer) { - checkFftBufferCounts(n: n, type: type, input: input, output: output, work: nil) - pffft_zreorder(ptr, input.baseAddress, output.baseAddress, PFFFT_FORWARD) - } - - public func convolveAccumulate(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) - } - - public func convolve(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) - } - - public static func minFftSize(for type: FFTType) -> Int { - Int(pffft_min_fft_size(pffft_transform_t(type))) - } - - public static func isValidSize(_ n: Int, for type: FFTType) -> Bool { - pffft_is_valid_size(Int32(n), pffft_transform_t(type)) != 0 - } - - 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)) - } - - deinit { - pffft_destroy_setup(ptr) - } -} - -extension Float: FFTElement { - public typealias FFTImpl = FFTFloat -} diff --git a/Packages/PFFFT/Sources/PFFFT/PFFFT.swift b/Packages/PFFFT/Sources/PFFFT/PFFFT.swift index 007f6d5..f2a5b1d 100644 --- a/Packages/PFFFT/Sources/PFFFT/PFFFT.swift +++ b/Packages/PFFFT/Sources/PFFFT/PFFFT.swift @@ -1,27 +1,12 @@ internal import PFFFTLib +import ComplexModule +import RealModule -@frozen -public enum FFTType { - case real - case complex -} -@frozen -public enum FFTSign: Int { - case forward = -1 - case backward = 1 -} - -public enum FFTError: Error { - case invalidSize -} - -public protocol FFTElement { - associatedtype FFTImpl: PFFFTProtocol -} - -public protocol PFFFTProtocol { - associatedtype Element +public protocol PFFFTProtocol { + associatedtype T: FFTElemental + associatedtype ComplexType = T.ComplexType + associatedtype ScalarType = T.ScalarType /// Get an FFT interface for the given size and FFT Type. /// @@ -31,14 +16,14 @@ public protocol PFFFTProtocol { /// - type: The type of FFT. /// - Returns: An FFT interface. /// - Throws: `FFTError.invalidSize` if the size is invalid. - static func setup(for n: Int, type: FFTType) throws -> Self + static func setup(for n: Int) throws -> Self /// 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 + init(n: Int) throws /// Perform a forward FFT on the input buffer. /// @@ -60,9 +45,9 @@ public protocol PFFFTProtocol { /// - 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 forward(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?) + func forward(input: borrowing Buffer, output: borrowing Buffer) - func inverse(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?) + func inverse(input: borrowing Buffer, output: borrowing Buffer) /// Perform a forward FFT on the input buffer, with implementation defined order. /// @@ -73,22 +58,11 @@ public protocol PFFFTProtocol { /// - 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 forwardUnordered(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?) + func forwardToInternalLayout(input: borrowing Buffer, output: borrowing Buffer) - func inverseUnordered(input: borrowing Buffer, output: borrowing Buffer, work: borrowing Buffer?) + func inverseFromInternalLayout(input: borrowing Buffer, output: borrowing Buffer) - func reorder(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 convolveAccumulate(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Element) + func reorder(spectrum: borrowing Buffer, output: borrowing Buffer) /// Perform a convolution of two complex signals in the frequency domain. /// @@ -99,12 +73,24 @@ public protocol PFFFTProtocol { /// - 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 convolve(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: Element) + func convolve(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: ScalarType) + + /// 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 convolveAccumulate(dftA: borrowing Buffer, dftB: borrowing Buffer, dftAB: borrowing Buffer, scaling: ScalarType) + /// Returns the minimum FFT size for the given type. /// /// - Parameter type: The type of FFT. - static func minFftSize(for type: FFTType) -> Int + static func minFftSize() -> Int /// Returns whether the given size is valid for the given type. /// @@ -112,10 +98,10 @@ public protocol PFFFTProtocol { /// - Parameters: /// - n: The size to check. /// - type: The type of FFT. - static func isValidSize(_ n: Int, for type: FFTType) -> Bool + static func isValidSize(_ n: Int) -> Bool /// Returns the nearest valid size for the given type. - static func nearestValidSize(_ n: Int, for type: FFTType, higher: Bool) -> Int + static func nearestValidSize(_ n: Int, higher: Bool) -> Int } public var simdArch: String { @@ -142,37 +128,3 @@ 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 10bbc29..8bb3c73 100644 --- a/Packages/PFFFT/Sources/PFFFT/SetupCache.swift +++ b/Packages/PFFFT/Sources/PFFFT/SetupCache.swift @@ -1,28 +1,24 @@ import Foundation -struct CacheKey : Hashable { - let n: Int - let type: FFTType -} - -public class SetupCache : @unchecked Sendable { - private var cache: [CacheKey: Element.FFTImpl?] = [:] +public class SetupCache : @unchecked Sendable { + typealias FFTType = FFT + private var cache: [Int: FFT] = [:] private let queue = DispatchQueue(label: String(describing: SetupCache.self), attributes: .concurrent) public init() {} - public func get(for n: Int, type: FFTType) throws -> Element.FFTImpl { - var setup: Element.FFTImpl?? + public func get(for n: Int) throws -> FFT { + var setup: FFTType?? queue.sync { - setup = cache[CacheKey(n: n, type: type)] + setup = cache[n] } if setup == nil { queue.sync(flags: .barrier) { - setup = cache[CacheKey(n: n, type: type)] + setup = cache[n] if setup == nil { - let entry = try? Element.FFTImpl(n: n, type: type) - cache[CacheKey(n: n, type: type)] = entry + let entry = try? FFTType(n: n) + cache[n] = entry setup = entry } } diff --git a/Sources/EcgSynKit/EcgSynKit.swift b/Sources/EcgSynKit/EcgSynKit.swift index 96e2860..a7cad99 100644 --- a/Sources/EcgSynKit/EcgSynKit.swift +++ b/Sources/EcgSynKit/EcgSynKit.swift @@ -1,7 +1,9 @@ import Algorithms +import ComplexModule import Foundation -import PFFFTLib import PFFFT +import PFFFTLib +import RealModule struct Parameters { /// The number of beats to simulate. @@ -63,75 +65,70 @@ func stdev(_ data: [Double]) -> Double { return sqrt(data.lazy.map { ($0 - mean) * ($0 - mean) }.reduce(0.0, +) / (n - 1)) } -struct RRProcess : ~Copyable { +struct RRProcess: ~Copyable { let nrr: Int - let input: Buffer - let output: Buffer - let work: Buffer - let fft: FFTDouble + let spectrum: Buffer> + let signal: Buffer + let fft: FFT init(nrr: Int) { self.nrr = nrr - self.input = Buffer(capacity: nrr + 2) - self.output = Buffer(capacity: nrr) - self.work = Buffer(capacity: nrr) - self.fft = try! FFTDouble.setup(for: nrr, type: .real) + fft = try! FFT(n: nrr) + spectrum = fft.makeSpectrumBuffer() + signal = fft.makeSignalBuffer() } func generate(params: Parameters) -> [Double] { + let w1 = 2.0 * .pi * params.flo + let w2 = 2.0 * .pi * params.fhi + let c1 = 2.0 * .pi * params.flostd + let c2 = 2.0 * .pi * params.fhistd - let w1 = 2.0 * .pi * params.flo - let w2 = 2.0 * .pi * params.fhi - let c1 = 2.0 * .pi * params.flostd - let c2 = 2.0 * .pi * params.fhistd + let sig2 = 1.0 + let sig1 = params.lfhfRatio - let sig2 = 1.0 - let sig1 = params.lfhfRatio + let rrmean = 60.0 / params.hrMean + let rrstd = 60.0 * params.hrStd / (params.hrMean * params.hrMean) - let rrmean = 60.0 / params.hrMean - let rrstd = 60.0 * params.hrStd / (params.hrMean * params.hrMean) + let sf = Double(params.sfInternal) + let df = sf / Double(nrr) - let sf = Double(params.sfInternal) - let df = sf / Double(nrr) + spectrum.withUnsafeMutableBufferPointer { swc in + for i in 0 ..< nrr / 2 + 1 { + let w = df * Double(i) * 2.0 * .pi + let dw1 = w - w1 + let dw2 = w - w2 + let hw = sig1 * exp(-dw1 * dw1 / (2.0 * c1 * c1)) / sqrt(2.0 * .pi * c1 * c1) + + sig2 * exp(-dw2 * dw2 / (2.0 * c2 * c2)) / sqrt(2.0 * .pi * c2 * c2) - input.withUnsafeMutableBufferPointer { swc in - for i in 0 ..< nrr / 2 + 1 { - let w = df * Double(i) * 2.0 * .pi - let dw1 = w - w1 - let dw2 = w - w2 - let hw = sig1 * exp(-dw1 * dw1 / (2.0 * c1 * c1)) / sqrt(2.0 * .pi * c1 * c1) - + sig2 * exp(-dw2 * dw2 / (2.0 * c2 * c2)) / sqrt(2.0 * .pi * c2 * c2) + let sw = (sf / 2.0) * sqrt(hw) + let ph = 2.0 * .pi * Double.random(in: 0.0 ..< 1.0) - let sw = (sf / 2.0) * sqrt(hw) - let ph = 2.0 * .pi * Double.random(in: 0.0 ..< 1.0) + swc[i].real = sw * cos(ph) + swc[i].imaginary = sw * sin(ph) + } - swc[i * 2] = sw * cos(ph) - swc[i * 2 + 1] = sw * sin(ph) + // pack Nyquist frequency real to imaginary of DC + swc[0].imaginary = swc[nrr / 2].real } - // pack Nyquist frequency real to imaginary of DC - swc[1] = swc[nrr] - } + fft.inverse(spectrum: spectrum, signal: signal) - self.fft.forward(input: input, output: output, work: nil, sign: .backward) + var rr = signal.withUnsafeMutableBufferPointer { outptr in + 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) - let ratio = rrstd / xstd - - for i in 0 ..< nrr { - rr[i] = rr[i] * ratio + rrmean - } - return rr + let xstd = stdev(rr) + let ratio = rrstd / xstd + for i in 0 ..< nrr { + rr[i] = rr[i] * ratio + rrmean + } + return rr } } - // func compute(params: Parameters) { // // adjust extrema parameters for mean heart rate // let hrFact = sqrt(params.hrMean / 60) diff --git a/Tests/EcgSynKitTests/EcgSynKitTests.swift b/Tests/EcgSynKitTests/EcgSynKitTests.swift index e18059b..2df535f 100644 --- a/Tests/EcgSynKitTests/EcgSynKitTests.swift +++ b/Tests/EcgSynKitTests/EcgSynKitTests.swift @@ -6,7 +6,8 @@ import Testing let p = Parameters() - let rr = rrprocess(params: p, nrr: 131072) + let rrp = RRProcess(nrr: 131072) + let rr = rrp.generate(params: p) //fft(data: &a, isign: -1) //print("out: \(rr)")