Swiftで高速フーリエ変換

slow_fft.swift

import Foundation

struct Complex64 {
    var real: Float
    var imag: Float

    init(_ r: Float, _ i: Float) {
        self.real = r
        self.imag = i
    }
}

extension Complex64 {
    static func +(left: Complex64, right: Complex64) -> Complex64 {
        return Complex64(left.real + right.real, left.imag + right.imag)
    }
    static func +=(left: inout Complex64, right: Complex64) {
        left = left + right
    }
    static func -(left: Complex64, right: Complex64) -> Complex64 {
        return Complex64(left.real - right.real, left.imag - right.imag)
    }
    static func -=(left: inout Complex64, right: Complex64) {
        left = left - right
    }
    static func *(left: Complex64, right: Complex64) -> Complex64 {
        return Complex64(left.real * right.real - left.imag * right.imag, left.real * right.imag + left.imag * right.real)
    }
    static func *=(left: inout Complex64, right: Complex64) {
        left = left * right
    }
}

func getIndex(_ length: Int) -> [Int] {
    let pow = Int(log2(Float(length)))
    var index: [Int] = []

    index.append(0)
    index.append(1)

    for _ in 0 ..< (pow - 1) {
        let indexLength = index.count

        for j in 0 ..< index.count {
            index[j] *= 2
        }
        for j in 0 ..< indexLength {
            index.append(index[j])
        }
        for j in indexLength ..< index.count {
            index[j] += 1
        }
    }

    return index
}

func FFT(signal: [Float]) -> [Complex64] {
    let length = signal.count
    let index = getIndex(length)
    let pow = Int(log2(Float(length)))

    var output = [Complex64]()

    for i in 0 ..< length {
        output.append(Complex64(signal[index[i]], 0.0))
    }

    var po2 = 1

    for _ in 1 ... pow {
        po2 = po2 << 1

        let po2m = po2 >> 1
        let theta = 2.0 * Double.pi / Double(po2)
        let w = Complex64(Float(cos(theta)), -Float(sin(theta)))

        var ws = Complex64(1.0, 0.0)

        for k in 0 ..< po2m {
            for j in stride(from: 0, to: length, by: po2) {
                let wfb = ws * output[j+k+po2m]

                output[j+k+po2m] = output[j+k] - wfb
                output[j+k] += wfb
            }

            ws *= w
        }
    }

    return output
}

let signal = [Float](repeating: 0.0, count: 32768)
let start = Date()
let _ = FFT(signal: signal)
let elapsed = Date().timeIntervalSince(start)

print(elapsed)

vdsp_fft.swift

import Accelerate

struct Complex64 {
    var real: Float
    var imag: Float

    init(_ r: Float, _ i: Float) {
        self.real = r
        self.imag = i
    }
}

func FFT(signal: [Float]) -> [Complex64] {
    let log2n = vDSP_Length(log2(Float(signal.count)) + 1)
    let fft = vDSP.FFT(log2n: log2n, radix: .radix2, ofType: DSPSplitComplex.self)!

    var inputReal = [Float](signal)
    var inputImag = [Float](repeating: 0.0, count: signal.count)
    var outputReal = [Float](repeating: 0.0, count: signal.count)
    var outputImag = [Float](repeating: 0.0, count: signal.count)

    inputReal.withUnsafeMutableBufferPointer { inputRealPtr in
        inputImag.withUnsafeMutableBufferPointer { inputImagPtr in
            outputReal.withUnsafeMutableBufferPointer { outputRealPtr in
                outputImag.withUnsafeMutableBufferPointer { outputImagPtr in
                    let input = DSPSplitComplex(realp: inputRealPtr.baseAddress!, imagp: inputImagPtr.baseAddress!)
                    var output = DSPSplitComplex(realp: outputRealPtr.baseAddress!, imagp: outputImagPtr.baseAddress!)

                    fft.forward(input: input, output: &output)
                }
            }
        }
    }

    var spectrum = [Complex64](repeating: Complex64(0.0, 0.0), count: signal.count)

    for i in 0 ..< signal.count {
        spectrum[i].real = outputReal[i]
        spectrum[i].imag = outputImag[i]
    }

    return spectrum
}

let signal = [Float](repeating: 0.0, count: 32768)
let start = Date()
let _ = FFT(signal: signal)
let elapsed = Date().timeIntervalSince(start)

print(elapsed)