はじめに

Accelerateフレームワークを利用すると数値計算やデジタル信号処理を高速かつ省電力で実行できます。vDSPはAccelerateフレームワークに含まれる機能の一つであり、デジタル信号処理を担当します。

vDSPには1次元の高速フーリエ変換を行うvDSP.FFTが含まれています。この機能はiOS 14から利用可能です。AccelerateフレームワークはAppleが提供する標準フレームワークであり追加のインストール作業などは不要です。

vDSPを利用したFFT

vDSPを利用したFFTの実装例を示します。

import Accelerate

// ダミーの信号を生成する関数です。周波数と振幅のペアが格納された配列を受け取り、複数のサイン波が合成された信号を生成します。
func synthesizeSignal(frequencyAmplitudePairs: [(f: Float, a: Float)], count: Int) -> [Float] {
  let tau: Float = .pi * 2
  let signal: [Float] = (0..<count).map { index in
    frequencyAmplitudePairs.reduce(0) { accumulator, frequenciesAmplitudePair in
      let normalizedIndex = Float(index) / Float(count)
      return accumulator + sin(normalizedIndex * frequenciesAmplitudePair.f * tau)
        * frequenciesAmplitudePair.a
    }
  }

  return signal
}

// FFTのサンプル数は2048とします。
let n = vDSP_Length(2048)

let frequencyAmplitudePairs = [
  (f: Float(2), a: Float(0.8)),
  (f: Float(7), a: Float(1.2)),
  (f: Float(24), a: Float(0.7)),
  (f: Float(50), a: Float(1.0)),
]

// 周波数と振幅のペアからダミーの信号を生成します。
let signal = synthesizeSignal(frequencyAmplitudePairs: frequencyAmplitudePairs, count: Int(n))

// FFTのサンプル数は指数で指示します。
let log2n = vDSP_Length(log2(Float(n)))

// 効率的にFFTの順方向・逆方向計算をするため、事前に三角関数テーブルを作成します。
guard let fftSetUp = vDSP.FFT(log2n: log2n, radix: .radix2, ofType: DSPSplitComplex.self) else {
  fatalError("Can't create FFT Setup.")
}

// 信号はナイキスト周波数で折り返されるため、入出力の配列の長さを`FFTのサンプル数 / 2`とします。
let halfN = Int(n / 2)

// 入力と出力について、複素数の実部と虚部を別々の配列として作成します。
var forwardInputReal = [Float](repeating: 0, count: halfN)
var forwardInputImag = [Float](repeating: 0, count: halfN)
var forwardOutputReal = [Float](repeating: 0, count: halfN)
var forwardOutputImag = [Float](repeating: 0, count: halfN)

forwardInputReal.withUnsafeMutableBufferPointer { forwardInputRealPtr in
  forwardInputImag.withUnsafeMutableBufferPointer { forwardInputImagPtr in
    forwardOutputReal.withUnsafeMutableBufferPointer { forwardOutputRealPtr in
      forwardOutputImag.withUnsafeMutableBufferPointer { forwardOutputImagPtr in
        // 入力信号を格納する`DSPSplitComplex`を作成します。
        var forwardInput = DSPSplitComplex(
          realp: forwardInputRealPtr.baseAddress!, imagp: forwardInputImagPtr.baseAddress!)

        // `signal`の実数値を複素数に変換します。
        signal.withUnsafeBytes {
          vDSP.convert(
            interleavedComplexVector: [DSPComplex]($0.bindMemory(to: DSPComplex.self)),
            toSplitComplexVector: &forwardInput)
        }

        // FFTの結果を受け取るために`DSPSplitComplex`を作成します。
        var forwardOutput = DSPSplitComplex(
          realp: forwardOutputRealPtr.baseAddress!, imagp: forwardOutputImagPtr.baseAddress!)

        // 順方向のFFTを実行します。
        fftSetUp.forward(input: forwardInput, output: &forwardOutput)
      }
    }
  }
}

let autospectrum = [Float](unsafeUninitializedCapacity: halfN) {
  autospectrumBuffer, initializedCount in
  // `vDSP_zaspec`関数は出力を蓄積します。スペクトルを計算する前に初期化されていない`autospectrumBuffer`をクリアします。
  vDSP.clear(&autospectrumBuffer)

  forwardOutputReal.withUnsafeMutableBufferPointer { forwardOutputRealPtr in
    forwardOutputImag.withUnsafeMutableBufferPointer { forwardOutputImagPtr in
      var frequencyDomain = DSPSplitComplex(
        realp: forwardOutputRealPtr.baseAddress!, imagp: forwardOutputImagPtr.baseAddress!)

      vDSP_zaspec(&frequencyDomain, autospectrumBuffer.baseAddress!, vDSP_Length(halfN))
    }
  }

  initializedCount = halfN
}

// スペクトラムから周波数と振幅のペアを求めます。
let componentFrequencyAmplitudePairs = autospectrum.enumerated().filter {
  $0.element > 1
}.map {
  return ($0.offset, sqrt($0.element) / Float(n))
}

for pair in componentFrequencyAmplitudePairs {
  print("frequency: \(pair.0), amplitude: \(String(format: "%.2f", pair.1))")
}

上記はApple Developerドキュメントに掲載されているSwiftコードを一部変更して注釈のコメントを加えたものです。元の実装については以下を参照ください。

Finding the component frequencies in a composite sine wave - Apple Developer Documentation

（おまけ）低速なFFT

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

高速フーリエ変換(FFT)をおじさんもC++で作ってみたよ - nursの日記

import Foundation

struct Complex64 {
  let real: Float
  let imag: Float
}

extension Complex64 {
  static func + (left: Complex64, right: Complex64) -> Complex64 {
    return Complex64(real: left.real + right.real, imag: left.imag + right.imag)
  }
  static func += (left: inout Complex64, right: Complex64) {
    left = left + right
  }
  static func - (left: Complex64, right: Complex64) -> Complex64 {
    return Complex64(real: left.real - right.real, imag: left.imag - right.imag)
  }
  static func -= (left: inout Complex64, right: Complex64) {
    left = left - right
  }
  static func * (left: Complex64, right: Complex64) -> Complex64 {
    return Complex64(
      real: left.real * right.real - left.imag * right.imag,
      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 forward(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(real: signal[index[i]], imag: 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(real: Float(cos(theta)), imag: -Float(sin(theta)))

    var ws = Complex64(real: 1.0, imag: 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
}

実行時間の比較

MacBook Pro 14インチ 2021年モデルで実行時間を計測しました。FFTのサンプルサイズは32768としました。以下の値はそれぞれ10回実行した中央値です。

結果

手作業で実装したFFT: 66.428 ms
vDSPを利用したFFT: 0.076 ms

およそ900倍の高速化が達成されています、やったね！

参考資料

Discussion