Add doc comments

This commit is contained in:
2024-10-28 03:29:52 -05:00
parent b6ed079d42
commit 4dda811a1c
8 changed files with 240 additions and 110 deletions

View File

@@ -29,7 +29,7 @@ let package = Package(
), ),
.testTarget( .testTarget(
name: "PFFFTTests", name: "PFFFTTests",
dependencies: ["PFFFTLib"] dependencies: ["PFFFT"]
), ),
] ]
) )

View File

@@ -3,10 +3,12 @@ let bufferAlignment = 32
@frozen @frozen
public struct Buffer<Element>: ~Copyable { public struct Buffer<Element>: ~Copyable {
let buffer: UnsafeMutableBufferPointer<Element> let buffer: UnsafeMutableBufferPointer<Element>
var count: Int { buffer.count }
var baseAddress: UnsafeMutablePointer<Element> { buffer.baseAddress! }
init(unsafeUninitializedCapacity: Int) { public init(capacity: Int) {
buffer = UnsafeMutableRawBufferPointer.allocate( buffer = UnsafeMutableRawBufferPointer.allocate(
byteCount: MemoryLayout<Element>.stride * unsafeUninitializedCapacity, byteCount: MemoryLayout<Element>.stride * capacity,
alignment: bufferAlignment alignment: bufferAlignment
).bindMemory(to: Element.self) ).bindMemory(to: Element.self)
} }
@@ -15,15 +17,15 @@ public struct Buffer<Element>: ~Copyable {
buffer.deallocate() buffer.deallocate()
} }
func withUnsafeMutableBufferPointer<R>(_ body: (UnsafeMutableBufferPointer<Element>) throws -> R) rethrows -> R { public func withUnsafeMutableBufferPointer<R>(_ body: (UnsafeMutableBufferPointer<Element>) throws -> R) rethrows -> R {
try body(buffer) try body(buffer)
} }
func withUnsafeBufferPointer<R>(_ body: (UnsafeBufferPointer<Element>) throws -> R) rethrows -> R { public func withUnsafeBufferPointer<R>(_ body: (UnsafeBufferPointer<Element>) throws -> R) rethrows -> R {
try body(UnsafeBufferPointer(buffer)) try body(UnsafeBufferPointer(buffer))
} }
func withUnsafeMutableBytes<R>(_ body: (UnsafeMutableRawBufferPointer) throws -> R) rethrows -> R { public func withUnsafeMutableBytes<R>(_ body: (UnsafeMutableRawBufferPointer) throws -> R) rethrows -> R {
try body(UnsafeMutableRawBufferPointer(buffer)) try body(UnsafeMutableRawBufferPointer(buffer))
} }
} }

View File

@@ -1,57 +1,69 @@
internal import PFFFTLib internal import PFFFTLib
struct FFTDoubleImpl: FFTImplProtocol { public final class FFTDoubleImpl: FFTImplProtocol {
let ptr: OpaquePointer let ptr: OpaquePointer
let n: Int
let type: FFTType
init?(n: Int, type: FFTType) { public init(n: Int, type: FFTType) throws {
guard let ptr = pffftd_new_setup(Int32(n), pffft_transform_t(type)) else { guard let ptr = pffftd_new_setup(Int32(n), pffft_transform_t(type)) else { throw FFTError.invalidSize }
return nil
}
self.ptr = ptr self.ptr = ptr
self.n = n
self.type = type
} }
func fft(input: borrowing Buffer<Double>, output: borrowing Buffer<Double>, work: borrowing Buffer<Double>?, sign: FFTSign) { public func fft(input: borrowing Buffer<Double>, output: borrowing Buffer<Double>, work: borrowing Buffer<Double>?, sign: FFTSign) {
let work: UnsafeMutablePointer<Double>! = switch work { checkFftBufferCounts(n: n, type: type, input: input, output: output, work: work)
case let .some(b): b.buffer.baseAddress
let workAddress: UnsafeMutablePointer<Double>! = switch work {
case let .some(b): b.baseAddress
case .none: nil case .none: nil
} }
pffftd_transform_ordered(ptr, input.baseAddress, output.baseAddress, workAddress, pffft_direction_t(sign))
pffftd_transform_ordered(ptr, input.buffer.baseAddress, output.buffer.baseAddress, work, pffft_direction_t(sign))
} }
func fftUnordered(input: borrowing Buffer<Double>, output: borrowing Buffer<Double>, work: borrowing Buffer<Double>?, sign: FFTSign) { public func fftUnordered(input: borrowing Buffer<Double>, output: borrowing Buffer<Double>, work: borrowing Buffer<Double>?, sign: FFTSign) {
let work: UnsafeMutablePointer<Double>! = switch work { checkFftBufferCounts(n: n, type: type, input: input, output: output, work: work)
case let .some(b): b.buffer.baseAddress
let workAddress: UnsafeMutablePointer<Double>! = switch work {
case let .some(b): b.baseAddress
case .none: nil 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<Double>, output: borrowing Buffer<Double>, sign: FFTSign) { public func zReorder(input: borrowing Buffer<Double>, output: borrowing Buffer<Double>, sign: FFTSign) {
pffftd_zreorder(ptr, input.buffer.baseAddress, output.buffer.baseAddress, pffft_direction_t(sign)) 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<Double>, dftB: borrowing Buffer<Double>, dftAB: borrowing Buffer<Double>, scaling: Double) { public func zConvolveAccumulate(dftA: borrowing Buffer<Double>, dftB: borrowing Buffer<Double>, dftAB: borrowing Buffer<Double>, scaling: Double) {
pffftd_zconvolve_accumulate(ptr, dftA.buffer.baseAddress, dftB.buffer.baseAddress, dftAB.buffer.baseAddress, scaling) 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<Double>, dftB: borrowing Buffer<Double>, dftAB: borrowing Buffer<Double>, scaling: Double) { public func zConvolve(dftA: borrowing Buffer<Double>, dftB: borrowing Buffer<Double>, dftAB: borrowing Buffer<Double>, scaling: Double) {
pffftd_zconvolve_no_accu(ptr, dftA.buffer.baseAddress, dftB.buffer.baseAddress, dftAB.buffer.baseAddress, scaling) 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))) 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 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)) Int(pffftd_nearest_transform_size(Int32(n), pffft_transform_t(type), higher ? 1 : 0))
} }
static func simdArch() -> String { deinit {
String(cString: pffftd_simd_arch()) pffftd_destroy_setup(ptr)
} }
} }
extension Double: FFTElement {
public typealias FFTImpl = FFTDoubleImpl
}

View File

@@ -1,57 +1,69 @@
internal import PFFFTLib internal import PFFFTLib
struct FFTFloatImpl: FFTImplProtocol { public final class FFTFloatImpl: FFTImplProtocol {
let ptr: OpaquePointer let ptr: OpaquePointer
let n: Int
let type: FFTType
init?(n: Int, type: FFTType) { public init(n: Int, type: FFTType) throws {
guard let ptr = pffft_new_setup(Int32(n), pffft_transform_t(type)) else { guard let ptr = pffft_new_setup(Int32(n), pffft_transform_t(type)) else { throw FFTError.invalidSize }
return nil
}
self.ptr = ptr self.ptr = ptr
self.n = n
self.type = type
} }
func fft(input: borrowing Buffer<Float>, output: borrowing Buffer<Float>, work: borrowing Buffer<Float>?, sign: FFTSign) { public func fft(input: borrowing Buffer<Float>, output: borrowing Buffer<Float>, work: borrowing Buffer<Float>?, sign: FFTSign) {
let work: UnsafeMutablePointer<Float>! = switch work { checkFftBufferCounts(n: n, type: type, input: input, output: output, work: work)
case let .some(b): b.buffer.baseAddress
let workAddress: UnsafeMutablePointer<Float>! = switch work {
case let .some(b): b.baseAddress
case .none: nil case .none: nil
} }
pffft_transform_ordered(ptr, input.baseAddress, output.baseAddress, workAddress, pffft_direction_t(sign))
pffft_transform_ordered(ptr, input.buffer.baseAddress, output.buffer.baseAddress, work, pffft_direction_t(sign))
} }
func fftUnordered(input: borrowing Buffer<Float>, output: borrowing Buffer<Float>, work: borrowing Buffer<Float>?, sign: FFTSign) { public func fftUnordered(input: borrowing Buffer<Float>, output: borrowing Buffer<Float>, work: borrowing Buffer<Float>?, sign: FFTSign) {
let work: UnsafeMutablePointer<Float>! = switch work { checkFftBufferCounts(n: n, type: type, input: input, output: output, work: work)
case let .some(b): b.buffer.baseAddress
let workAddress: UnsafeMutablePointer<Float>! = switch work {
case let .some(b): b.baseAddress
case .none: nil 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<Float>, output: borrowing Buffer<Float>, sign: FFTSign) { public func zReorder(input: borrowing Buffer<Float>, output: borrowing Buffer<Float>, sign: FFTSign) {
pffft_zreorder(ptr, input.buffer.baseAddress, output.buffer.baseAddress, pffft_direction_t(sign)) 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<Float>, dftB: borrowing Buffer<Float>, dftAB: borrowing Buffer<Float>, scaling: Float) { public func zConvolveAccumulate(dftA: borrowing Buffer<Float>, dftB: borrowing Buffer<Float>, dftAB: borrowing Buffer<Float>, scaling: Float) {
pffft_zconvolve_accumulate(ptr, dftA.buffer.baseAddress, dftB.buffer.baseAddress, dftAB.buffer.baseAddress, scaling) 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<Float>, dftB: borrowing Buffer<Float>, dftAB: borrowing Buffer<Float>, scaling: Float) { public func zConvolve(dftA: borrowing Buffer<Float>, dftB: borrowing Buffer<Float>, dftAB: borrowing Buffer<Float>, scaling: Float) {
pffft_zconvolve_no_accu(ptr, dftA.buffer.baseAddress, dftB.buffer.baseAddress, dftAB.buffer.baseAddress, scaling) 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))) 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 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)) Int(pffft_nearest_transform_size(Int32(n), pffft_transform_t(type), higher ? 1 : 0))
} }
static func simdArch() -> String { deinit {
String(cString: pffft_simd_arch()) pffft_destroy_setup(ptr)
} }
} }
extension Float: FFTElement {
public typealias FFTImpl = FFTFloatImpl
}

View File

@@ -1,29 +1,111 @@
internal import PFFFTLib internal import PFFFTLib
enum FFTType { @frozen
public enum FFTType {
case real case real
case complex case complex
} }
enum FFTSign: Int { @frozen
public enum FFTSign: Int {
case forward = -1 case forward = -1
case backward = 1 case backward = 1
} }
protocol FFTImplProtocol { public enum FFTError: Error {
case invalidSize
}
public protocol FFTElement {
associatedtype FFTImpl: FFTImplProtocol
}
public protocol FFTImplProtocol<Element>: ~Copyable {
associatedtype Element 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<Element>, output: borrowing Buffer<Element>, work: borrowing Buffer<Element>?, sign: FFTSign) func fft(input: borrowing Buffer<Element>, output: borrowing Buffer<Element>, work: borrowing Buffer<Element>?, 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<Element>, output: borrowing Buffer<Element>, work: borrowing Buffer<Element>?, sign: FFTSign) func fftUnordered(input: borrowing Buffer<Element>, output: borrowing Buffer<Element>, work: borrowing Buffer<Element>?, sign: FFTSign)
func zReorder(input: borrowing Buffer<Element>, output: borrowing Buffer<Element>, sign: FFTSign) func zReorder(input: borrowing Buffer<Element>, output: borrowing Buffer<Element>, 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<Element>, dftB: borrowing Buffer<Element>, dftAB: borrowing Buffer<Element>, scaling: Element) func zConvolveAccumulate(dftA: borrowing Buffer<Element>, dftB: borrowing Buffer<Element>, dftAB: borrowing Buffer<Element>, 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<Element>, dftB: borrowing Buffer<Element>, dftAB: borrowing Buffer<Element>, scaling: Element) func zConvolve(dftA: borrowing Buffer<Element>, dftB: borrowing Buffer<Element>, dftAB: borrowing Buffer<Element>, 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 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 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 { extension pffft_transform_t {
@@ -45,3 +127,38 @@ extension pffft_direction_t {
} }
} }
} }
@inline(__always)
func checkFftBufferCounts<Element>(n: Int, type: FFTType, input: borrowing Buffer<Element>, output: borrowing Buffer<Element>, work: borrowing Buffer<Element>?) {
let minSize = type == .real ? n : 2 * n
let work: UnsafeMutableBufferPointer<Element>! = 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<Element>(n: Int, type: FFTType, a: borrowing Buffer<Element>, b: borrowing Buffer<Element>, ab: borrowing Buffer<Element>) {
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")
}
}

View File

@@ -5,13 +5,15 @@ struct CacheKey : Hashable {
let type: FFTType let type: FFTType
} }
class SetupCache<Impl: FFTImplProtocol> : @unchecked Sendable { public class SetupCache<Element: FFTElement> : @unchecked Sendable {
private var cache: [CacheKey: Impl?] = [:] private var cache: [CacheKey: Element.FFTImpl?] = [:]
private let queue = DispatchQueue(label: String(describing: SetupCache.self), attributes: .concurrent) private let queue = DispatchQueue(label: String(describing: SetupCache.self), attributes: .concurrent)
func get(for n: Int, type: FFTType) -> Impl? { public init() {}
var setup: Impl??
public func get(for n: Int, type: FFTType) throws -> Element.FFTImpl {
var setup: Element.FFTImpl??
queue.sync { queue.sync {
setup = cache[CacheKey(n: n, type: type)] setup = cache[CacheKey(n: n, type: type)]
} }
@@ -19,12 +21,19 @@ class SetupCache<Impl: FFTImplProtocol> : @unchecked Sendable {
queue.sync(flags: .barrier) { queue.sync(flags: .barrier) {
setup = cache[CacheKey(n: n, type: type)] setup = cache[CacheKey(n: n, type: type)]
if setup == nil { if setup == nil {
setup = Impl(n: n, type: type) let entry = try? Element.FFTImpl(n: n, type: type)
cache[CacheKey(n: n, type: type)] = setup 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()
}
}
} }

View File

@@ -2,5 +2,14 @@ import Testing
@testable import PFFFT @testable import PFFFT
@Test func example() async throws { @Test func example() async throws {
let cache = PFFFT.SetupCache<Double>()
let ff = try cache.get(for: 1024, type: .real)
let input = Buffer<Double>(capacity: 1023)
let output = Buffer<Double>(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. // Write your test here and use APIs like `#expect(...)` to check expected conditions.
} }

View File

@@ -1,6 +1,7 @@
import Algorithms import Algorithms
import Foundation import Foundation
import PFFFTLib import PFFFTLib
import PFFFT
struct Parameters { struct Parameters {
/// The number of beats to simulate. /// The number of beats to simulate.
@@ -52,39 +53,6 @@ struct Parameters {
let fhistd = 0.01 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 { struct EcgDerive {
let rrpc: [Double] 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)) return sqrt(data.lazy.map { ($0 - mean) * ($0 - mean) }.reduce(0.0, +) / (n - 1))
} }
let fftlib = SetupCache<Double>()
func rrprocess(params: Parameters, nrr: Int) -> [Double] { func rrprocess(params: Parameters, nrr: Int) -> [Double] {
let w1 = 2.0 * .pi * params.flo let w1 = 2.0 * .pi * params.flo
let w2 = 2.0 * .pi * params.fhi let w2 = 2.0 * .pi * params.fhi
@@ -110,11 +80,11 @@ func rrprocess(params: Parameters, nrr: Int) -> [Double] {
let sf = Double(params.sfInternal) let sf = Double(params.sfInternal)
let df = sf / Double(nrr) let df = sf / Double(nrr)
let fft = PfftSetupCache.shared.get(for: nrr)
var rr = withUnsafeTemporaryAllocation(byteCount: (nrr + 2) * MemoryLayout<Double>.stride, alignment: 32) { let input = Buffer<Double>(capacity: (nrr + 2))
let swc = $0.bindMemory(to: Double.self) let output = Buffer<Double>(capacity: nrr)
input.withUnsafeMutableBufferPointer { swc in
for i in 0 ..< nrr / 2 + 1 { for i in 0 ..< nrr / 2 + 1 {
let w = df * Double(i) * 2.0 * .pi let w = df * Double(i) * 2.0 * .pi
let dw1 = w - w1 let dw1 = w - w1
@@ -131,13 +101,12 @@ func rrprocess(params: Parameters, nrr: Int) -> [Double] {
// pack Nyquist frequency real to imaginary of DC // pack Nyquist frequency real to imaginary of DC
swc[1] = swc[nrr] swc[1] = swc[nrr]
}
return withUnsafeTemporaryAllocation(byteCount: nrr * MemoryLayout<Double>.stride, alignment: 32) { try! fftlib.get(for: nrr, type: .real).fft(input: input, output: output, work: nil, sign: .backward)
let outptr = $0.bindMemory(to: Double.self)
pffftd_transform_ordered(fft, swc.baseAddress, outptr.baseAddress, nil, PFFFT_BACKWARD) var rr = output.withUnsafeMutableBufferPointer { outptr in
return outptr.map { $0 * 1.0 / Double(nrr) } return outptr.map { $0 * 1.0 / Double(nrr) }
}
} }
let xstd = stdev(rr) let xstd = stdev(rr)