Zenn
😺

（Swift）vDSPを利用して高速フーリエ変換する

2021/03/13に公開
iOS
Swift
Apple
#
音声信号処理
#
vdsp
tech

vDSPとは

数値計算やデジタル信号処理を高速に行うためのライブラリがvDSPです。Xcodeに最初から含まれているため、追加インストールする必要はありません。

さて、vDSPには1次元の高速フーリエ変換を行うvDSP.FFTがあります。iOS 14から利用できるようになったので試してみましょう。

vDSPを利用した高速フーリエ変換

例として、複素数の実部のみが記録されたFloatの配列を受け取り、結果を複素数として返す処理を示します。音声信号などの解析を想定しています。

なお、Swiftは複素数型が提供されていないため独自にComplex64を定義しています。関数の戻り値として実部の配列と虚部の配列を返しても良いのですが、後でvDSPを使わない場合との比較をするため、このように実装しています。

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
}

低速な高速フーリエ変換

vDSPがどれほど高速なのか確かめるため、ライブラリを一切使わずに高速フーリエ変換を実装します。実装方法については以下の記事を参考にしました。

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
}

実行時間の比較

大雑把ではありますが、32768サンプルのデータをフーリエ変換した場合の実行時間を比較します。以下のように計測します。

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

print(elapsed)

計測に使用したSwiftのコードはgistで公開しています。

実はvDSP.FFTはmacOSにも提供されているため、全く同じコードがmacOSとiOSどちらでも動作します。実機のiOS 14端末を準備するのが面倒だったので、今回はswiftコマンドを使って計測することにします。

$ swift --version
Apple Swift version 5.3.2 (swiftlang-1200.0.45 clang-1200.0.32.28)
Target: x86_64-apple-darwin19.6.0

まずは自力で実装した低速な高速フーリエ変換から。

$ swift slow_fft.swift
0.13674700260162354

次はvDSPを利用した高速フーリエ変換です。

$ swift vdsp_fft.swift
0.022410988807678223

すばらしい、実行時間は137 msから22 msに短縮されました！約6倍の高速化です！

  • 注1 ... 実行時間の計測はそれぞれ10回実行して中央値を選びました。
  • 注2 ... 計測結果はMacBook Air 2020 Intel Core i5モデルのものです。

（補足）vDSPでフーリエ変換した場合のゼロ埋めについて

vDSPでフーリエ変換した結果は入力信号の配列の長さ / 2以降のデータがゼロ埋めされています。

配列の半分はエイリアシングで折り返された値です。通常は使われない値ですのでゼロ埋めされていても問題ないはずですが、その点は注意してください。

Discussion