Swiftで高速フーリエ変換（FFT）を実装する方法

iOSデバイスで高速フーリエ変換を実装する方法を紹介します。この記事では、まずはじめに一般的な高速フーリエ変換（FFT）であるCooley-Tukey法をSwiftで検証し、その後Accelerateというライブラリによる実装方法も併せて紹介します。

こんにちは。wat（@watlablog）です。この記事ではSwiftによる高速フーリエ変換について色々まとめてみます！　

目次（項目クリックでジャンプできます）

この記事の目標

Cooley-Tukey法によるFFTコードを自分で書いてみる

ここしばらく筆者はiPhoneアプリ制作をしてみたいと思い、基礎文法のチェックや録音の仕方、音声再生方法や画面遷移の方法を個別の記事で取り上げました。信号処理アプリを作る目標を立てており、この分野ではかかせない周波数分析機能に着手するのがこの記事の内容です。

WATLABブログでは以下の記事で前回離散フーリエ変換（DFT）をPythonやSwiftで実装してみました。
・Pythonで検証しながらSwiftで離散フーリエ変換を実装する

AccelerateでFFTを使ってみる

ただ、Pythonで実感したのは自分で実装したFFT関数よりもNumPy等の外部ライブラリを使った方が圧倒的に高速だということでした。基本的なアルゴリズムを学ぶために自分でコードを書いてみるのは非常に良いことですが、実用には既存のライブラリを使った方が良いと思います。

そのためこの記事では一度Cooley-Tukey法を書いてみた後に、SwiftのAccelerateというライブラリを使ってFFTを使う方法も学びます。

動作環境

この記事のコードは以下の開発環境で動作を確認しました。

Mac	OS	macOS Ventura 13.2.1
	CPU	1.4[GHz]
	メモリ	8[GB]
Xcode	Version	14.3(14E222b)
Swift	Version	5.8(swiftlang-5.8.0.124.2 clang-1403.0.22.11.100)
iPhone SE2	OS	iOS 16.3.1

SwiftでCooley-Tukey法によるFFTをフルスクラッチするコード

事前準備

まずは自作のFFTコードから書いてみますが、動作検証のためにいくつか事前準備を行います。

データ点数を指定して配列を作成する関数

FFTのアルゴリズムを効果的に使用するためにはデータ点数が2の冪乗でないとなりません。Swiftでメソッドがあるかどうか探せなかったので、まずはデータ点数を自由に指定するために以下のlinspace関数を作りました。Pythonのnp.linspace()を真似たのですが、これでFFTコードの検証がしやすくなります。

func linspace(start: Float, stop: Float, num: Int) -> [Float] {
        // データ点数を指定して配列を作成する関数
        
        let increment = (stop - start) / Float(num - 1)
        return (0..<num).map { Float($0) * increment + start }
    }

func linspace(start: Float, stop: Float, num: Int) -> [Float] {

// データ点数を指定して配列を作成する関数

let increment = (stop - start) / Float(num - 1)

return (0..<num).map { Float($0) * increment + start }

}

サンプル波形を生成する関数

上記データ点数を指定するlinspace関数を使って、サイン波を作るwavePoints関数も用意しておきます。サイン波は周波数や時間長を変更できます。また、コードの検証としては直流成分も確認したいため、dcで0[Hz]の時の振幅を定義できるようにしています。

func wavePoints(N: Int, length: Float) -> (xArray: [Float], yArray: [Float]){
        // プロットするデータを演算する関数
        
        // データ点数Nと時間長lengthからx軸を作成
        let xArray = linspace(start: 0.0, stop: length, num: N)
        
        // 正弦波を作成。A:振幅、f:周波数[Hz]、dc:直流成分
        let A: Float = 1
        let f: Float = 100
        let dc: Float = 1
        let yArray = xArray.map{ A * sin(2 * .pi * f * $0) + dc }
        
        return (xArray, yArray)
    }

func wavePoints(N: Int, length: Float) -> (xArray: [Float], yArray: [Float]){

// プロットするデータを演算する関数

// データ点数Nと時間長lengthからx軸を作成

let xArray = linspace(start: 0.0, stop: length, num: N)

// 正弦波を作成。A:振幅、f:周波数[Hz]、dc:直流成分

let A: Float = 1

let f: Float = 100

let dc: Float = 1

let yArray = xArray.map{ A * sin(2 * .pi * f * $0) + dc }

return (xArray, yArray)

}

複素数構造体を定義

Swiftには標準で複素数型がありません。そのためstructを使ってComplex構造体を定義します。詳細は離散フーリエ変換の記事で書いた通りです。

struct Complex {
        // 複素数型の構造体を定義
        
        // re:実部、im:虚部
        var re: Float
        var im: Float
        
        // Complex + Complexの場合の演算処理
        static func +(lhs: Complex, rhs: Complex) -> Complex {
            return Complex(re: lhs.re + rhs.re, im: lhs.im + rhs.im)
        }
        
        // Complex * Complexの場合の演算処理
        static func *(lhs: Complex, rhs: Complex) -> Complex {
            return Complex(re: lhs.re * rhs.re - lhs.im * rhs.im, im: lhs.re * rhs.im + lhs.im * rhs.re)
        }
        
        // Complex * Floatの場合の演算処理
        static func *(lhs: Complex, rhs: Float) -> Complex {
            return Complex(re: lhs.re * rhs, im: lhs.im * rhs)
        }
    }

struct Complex {

// 複素数型の構造体を定義

// re:実部、im:虚部

var re: Float

var im: Float

// Complex + Complexの場合の演算処理

static func +(lhs: Complex, rhs: Complex) -> Complex {

return Complex(re: lhs.re + rhs.re, im: lhs.im + rhs.im)

}

// Complex * Complexの場合の演算処理

static func *(lhs: Complex, rhs: Complex) -> Complex {

return Complex(re: lhs.re * rhs.re - lhs.im * rhs.im, im: lhs.re * rhs.im + lhs.im * rhs.re)

}

// Complex * Floatの場合の演算処理

static func *(lhs: Complex, rhs: Float) -> Complex {

return Complex(re: lhs.re * rhs, im: lhs.im * rhs)

}

FFT関数（Cooley-Tukey法）

こちらがSwiftで書いたCooley-Tukey法によるFFT関数です。先ほど紹介したPythonによるFFTの記事と同じ書き方をしているので、Pythonユーザーは是非両コードを見比べてみてください。

func fft(_ x: [Float]) -> [Complex] {
        // FFT Cooley-Tukey法
        
        // 標本点数
        let N = x.count
        
        // Nが1なら、そのままxを返す
        if N == 1 {
            return [Complex(re: x[0], im: 0)]
        }
        
        // ビット反転順序に基づいて、偶数インデックスと奇数インデックスの部分配列を得る
        var evenArray: [Float] = []
        var oddArray: [Float] = []
        for i in 0..<N {
            if i % 2 == 0 {
                evenArray.append(x[i])
            } else {
                oddArray.append(x[i])
            }
        }
        
        let even = fft(evenArray)
        let odd = fft(oddArray)
        
        // DFTの結果を格納するための複素数配列
        var T: [Complex] = []
        for k in 0..<N {
            let theta = -2.0 * Float.pi * Float(k) / Float(N)
            T.append(Complex(re: cos(theta), im: sin(theta)))
        }
        
        // 偶数部分と奇数部分のDFTを組み合わせる（バタフライ演算）
        var F: [Complex] = []
        for k in 0..<N {
            if k < N / 2 {
                F.append(even[k] + T[k] * odd[k])
            } else {
                F.append(even[k - N / 2] + T[k] * odd[k - N / 2])
            }
        }
        
        return F
    }

func fft(_ x: [Float]) -> [Complex] {

// FFT Cooley-Tukey法

// 標本点数

let N = x.count

// Nが1なら、そのままxを返す

if N == 1 {

return [Complex(re: x[0], im: 0)]

}

// ビット反転順序に基づいて、偶数インデックスと奇数インデックスの部分配列を得る

var evenArray: [Float] = []

var oddArray: [Float] = []

for i in 0..<N {

if i % 2 == 0 {

evenArray.append(x[i])

} else {

oddArray.append(x[i])

}

let even = fft(evenArray)

let odd = fft(oddArray)

// DFTの結果を格納するための複素数配列

var T: [Complex] = []

for k in 0..<N {

let theta = -2.0 * Float.pi * Float(k) / Float(N)

T.append(Complex(re: cos(theta), im: sin(theta)))

}

// 偶数部分と奇数部分のDFTを組み合わせる（バタフライ演算）

var F: [Complex] = []

for k in 0..<N {

if k < N / 2 {

F.append(even[k] + T[k] * odd[k])

} else {

F.append(even[k - N / 2] + T[k] * odd[k - N / 2])

}

return F

}

全コード（コピペ動作OK）

UI部分を含めた全コードを以下に示します。freqLimitはなくても動作しますが、グラフに表示させる周波数を絞ることを目的としています。グラフへの表示は以下の記事で紹介したChartを使っています。
・SwiftUI：iOS16から追加されたChartsでグラフを作成

import SwiftUI
import Charts

struct PointsData: Identifiable {
    // 点群データの構造体
    var xValue: Float
    var yValue: Float
    var id = UUID()
}

struct ContentView: View {
    
    // データを定義
    @State private var data: [PointsData] = []
    @State private var x: [Float] = []
    @State private var y: [Float] = []
    
    func linspace(start: Float, stop: Float, num: Int) -> [Float] {
        // データ点数を指定して配列を作成する関数
        
        let increment = (stop - start) / Float(num - 1)
        return (0..<num).map { Float($0) * increment + start }
    }
    
    func wavePoints(N: Int, length: Float) -> (xArray: [Float], yArray: [Float]){
        // プロットするデータを演算する関数
        
        // データ点数Nと時間長lengthからx軸を作成
        let xArray = linspace(start: 0.0, stop: length, num: N)
        
        // 正弦波を作成。A:振幅、f:周波数[Hz]、dc:直流成分
        let A: Float = 1
        let f: Float = 100
        let dc: Float = 1
        let yArray = xArray.map{ A * sin(2 * .pi * f * $0) + dc }
        
        return (xArray, yArray)
    }
    
    struct Complex {
        // 複素数型の構造体を定義
        
        // re:実部、im:虚部
        var re: Float
        var im: Float
        
        // Complex + Complexの場合の演算処理
        static func +(lhs: Complex, rhs: Complex) -> Complex {
            return Complex(re: lhs.re + rhs.re, im: lhs.im + rhs.im)
        }
        
        // Complex * Complexの場合の演算処理
        static func *(lhs: Complex, rhs: Complex) -> Complex {
            return Complex(re: lhs.re * rhs.re - lhs.im * rhs.im, im: lhs.re * rhs.im + lhs.im * rhs.re)
        }
        
        // Complex * Floatの場合の演算処理
        static func *(lhs: Complex, rhs: Float) -> Complex {
            return Complex(re: lhs.re * rhs, im: lhs.im * rhs)
        }
    }
    
    func fft(_ x: [Float]) -> [Complex] {
        // FFT Cooley-Tukey法
        
        // 標本点数
        let N = x.count
        
        // Nが1なら、そのままxを返す
        if N == 1 {
            return [Complex(re: x[0], im: 0)]
        }
        
        // ビット反転順序に基づいて、偶数インデックスと奇数インデックスの部分配列を得る
        var evenArray: [Float] = []
        var oddArray: [Float] = []
        for i in 0..<N {
            if i % 2 == 0 {
                evenArray.append(x[i])
            } else {
                oddArray.append(x[i])
            }
        }
        
        let even = fft(evenArray)
        let odd = fft(oddArray)
        
        // DFTの結果を格納するための複素数配列
        var T: [Complex] = []
        for k in 0..<N {
            let theta = -2.0 * Float.pi * Float(k) / Float(N)
            T.append(Complex(re: cos(theta), im: sin(theta)))
        }
        
        // 偶数部分と奇数部分のDFTを組み合わせる（バタフライ演算）
        var F: [Complex] = []
        for k in 0..<N {
            if k < N / 2 {
                F.append(even[k] + T[k] * odd[k])
            } else {
                F.append(even[k - N / 2] + T[k] * odd[k - N / 2])
            }
        }
        
        return F
    }
    
    var body: some View {
        // UI
        
        VStack{
            Chart {
                // データ構造からx, y値を取得してプロット
                
                ForEach(data) { shape in
                    LineMark(
                        x: .value("x", shape.xValue),
                        y: .value("y", shape.yValue)
                    )
                }
            }
            .padding(.all, 10)
            
            Button("Button"){
                // ボタンイベント
                
                // 関数を実行してx, y演算結果を取得
                let result = wavePoints(N: 8192, length: 10)
                let x = result.xArray
                let y = result.yArray
                
                // FFTを計算
                let fftResult = fft(y)
                
                // FFTの結果から振幅を計算
                var amp: [Float] = []
                for i in 0..<fftResult.count {
                    let re = fftResult[i].re
                    let im = fftResult[i].im
                    let magnitude = sqrt(re * re + im * im)
                    if i == 0 {
                        // 直流成分
                        amp.append(magnitude / Float(fftResult.count))
                    } else {
                        // 変動成分
                        amp.append(2 * magnitude / Float(fftResult.count))
                    }
                }
                
                // 周波数軸を計算
                let dt = x[1] - x[0] // 時間刻み
                var freq: [Float] = []
                for k in 0..<x.count {
                    freq.append(Float(k) / (Float(x.count) * dt))
                }
                
                // 周波数の上限を指定
                let freqLimit: Float = 300.0
                let limitIndex = Int(freqLimit * Float(fftResult.count) * dt)
                
                let filteredFreq = freq[..<limitIndex]
                let filteredAmp = amp[..<limitIndex]
                
                // プロットデータの全削除
                data.removeAll()
                
                // プロットデータの追加
                for i in 0..<filteredFreq.count {
                    data.append(.init(xValue: filteredFreq[i], yValue: filteredAmp[i]))
                }
            }
        }
    }
}

struct ContentView_Previews: PreviewProvider {
    static var previews: some View {
        // プレビュー
        ContentView()
    }
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

import SwiftUI

import Charts

struct PointsData: Identifiable {

// 点群データの構造体

var xValue: Float

var yValue: Float

var id = UUID()

}

struct ContentView: View {

// データを定義

@State private var data: [PointsData] = []

@State private var x: [Float] = []

@State private var y: [Float] = []

func linspace(start: Float, stop: Float, num: Int) -> [Float] {

// データ点数を指定して配列を作成する関数

let increment = (stop - start) / Float(num - 1)

return (0..<num).map { Float($0) * increment + start }

}

func wavePoints(N: Int, length: Float) -> (xArray: [Float], yArray: [Float]){

// プロットするデータを演算する関数

// データ点数Nと時間長lengthからx軸を作成

let xArray = linspace(start: 0.0, stop: length, num: N)

// 正弦波を作成。A:振幅、f:周波数[Hz]、dc:直流成分

let A: Float = 1

let f: Float = 100

let dc: Float = 1

let yArray = xArray.map{ A * sin(2 * .pi * f * $0) + dc }

return (xArray, yArray)

}

struct Complex {

// 複素数型の構造体を定義

// re:実部、im:虚部

var re: Float

var im: Float

// Complex + Complexの場合の演算処理

static func +(lhs: Complex, rhs: Complex) -> Complex {

return Complex(re: lhs.re + rhs.re, im: lhs.im + rhs.im)

}

// Complex * Complexの場合の演算処理

static func *(lhs: Complex, rhs: Complex) -> Complex {

return Complex(re: lhs.re * rhs.re - lhs.im * rhs.im, im: lhs.re * rhs.im + lhs.im * rhs.re)

}

// Complex * Floatの場合の演算処理

static func *(lhs: Complex, rhs: Float) -> Complex {

return Complex(re: lhs.re * rhs, im: lhs.im * rhs)

}

func fft(_ x: [Float]) -> [Complex] {

// FFT Cooley-Tukey法

// 標本点数

let N = x.count

// Nが1なら、そのままxを返す

if N == 1 {

return [Complex(re: x[0], im: 0)]

}

// ビット反転順序に基づいて、偶数インデックスと奇数インデックスの部分配列を得る

var evenArray: [Float] = []

var oddArray: [Float] = []

for i in 0..<N {

if i % 2 == 0 {

evenArray.append(x[i])

} else {

oddArray.append(x[i])

}

let even = fft(evenArray)

let odd = fft(oddArray)

// DFTの結果を格納するための複素数配列

var T: [Complex] = []

for k in 0..<N {

let theta = -2.0 * Float.pi * Float(k) / Float(N)

T.append(Complex(re: cos(theta), im: sin(theta)))

}

// 偶数部分と奇数部分のDFTを組み合わせる（バタフライ演算）

var F: [Complex] = []

for k in 0..<N {

if k < N / 2 {

F.append(even[k] + T[k] * odd[k])

} else {

F.append(even[k - N / 2] + T[k] * odd[k - N / 2])

}

return F

}

var body: some View {

// UI

VStack{

Chart {

// データ構造からx, y値を取得してプロット

ForEach(data) { shape in

LineMark(

x: .value("x", shape.xValue),

y: .value("y", shape.yValue)

)

}

.padding(.all, 10)

Button("Button"){

// ボタンイベント

// 関数を実行してx, y演算結果を取得

let result = wavePoints(N: 8192, length: 10)

let x = result.xArray

let y = result.yArray

// FFTを計算

let fftResult = fft(y)

// FFTの結果から振幅を計算

var amp: [Float] = []

for i in 0..<fftResult.count {

let re = fftResult[i].re

let im = fftResult[i].im

let magnitude = sqrt(re * re + im * im)

if i == 0 {

// 直流成分

amp.append(magnitude / Float(fftResult.count))

} else {

// 変動成分

amp.append(2 * magnitude / Float(fftResult.count))

}

// 周波数軸を計算

let dt = x[1] - x[0] // 時間刻み

var freq: [Float] = []

for k in 0..<x.count {

freq.append(Float(k) / (Float(x.count) * dt))

}

// 周波数の上限を指定

let freqLimit: Float = 300.0

let limitIndex = Int(freqLimit * Float(fftResult.count) * dt)

let filteredFreq = freq[..<limitIndex]

let filteredAmp = amp[..<limitIndex]

// プロットデータの全削除

data.removeAll()

// プロットデータの追加

for i in 0..<filteredFreq.count {

data.append(.init(xValue: filteredFreq[i], yValue: filteredAmp[i]))

}

struct ContentView_Previews: PreviewProvider {

static var previews: some View {

// プレビュー

ContentView()

}

実行結果

下図が実行結果の例です。Buttonをタップすると波形が表示されます。これはwavePointsで生成された時間波形をfft関数で高速フーリエ変換した結果です。指定した通りのピークがでました。

freqLimitを外したコード、もしくはfreqLimitの上限を最大まで持っていった場合は下図になります。fft関数で計算された値は鏡像ピークまで入っている結果です。

SwiftのAccelerate/vDSPでFFTをするコード

全コード

それではAccelerateというライブラリのvDSPを使ってみます。上記のコードをvDSPを使って書いたコードを以下に示します。

import SwiftUI
import Charts
import Accelerate

struct PointsData: Identifiable {
    // 点群データの構造体
    var xValue: Float
    var yValue: Float
    var id = UUID()
}

struct ContentView: View {
    
    // データを定義
    @State private var data: [PointsData] = []
    @State private var x: [Float] = []
    @State private var y: [Float] = []
    @State private var showingAlert = false
    
    func linspace(start: Float, stop: Float, num: Int) -> [Float] {
        // データ点数を指定して配列を作成する関数
        
        let increment = (stop - start) / Float(num - 1)
        return (0..<num).map { Float($0) * increment + start }
    }
    
    func wavePoints(N: Int, length: Float) -> (xArray: [Float], yArray: [Float]){
        // プロットするデータを演算する関数
        
        // データ点数Nと時間長lengthからx軸を作成
        let xArray = linspace(start: 0.0, stop: length, num: N)
        
        // 正弦波を作成。A:振幅、f:周波数[Hz]、dc:直流成分
        let A: Float = 1
        let f: Float = 100
        let dc: Float = 1
        let yArray = xArray.map{ A * sin(2 * .pi * f * $0) + dc }
        
        return (xArray, yArray)
    }
    
    struct Complex {
        // 複素数型の構造体を定義
        
        // re:実部、im:虚部
        var re: Float
        var im: Float
        
        // Complex + Complexの場合の演算処理
        static func +(lhs: Complex, rhs: Complex) -> Complex {
            return Complex(re: lhs.re + rhs.re, im: lhs.im + rhs.im)
        }
        
        // Complex * Complexの場合の演算処理
        static func *(lhs: Complex, rhs: Complex) -> Complex {
            return Complex(re: lhs.re * rhs.re - lhs.im * rhs.im, im: lhs.re * rhs.im + lhs.im * rhs.re)
        }
        
        // Complex * Floatの場合の演算処理
        static func *(lhs: Complex, rhs: Float) -> Complex {
            return Complex(re: lhs.re * rhs, im: lhs.im * rhs)
        }
    }
    
    func isPowerOfTwo(_ n: Int) -> Bool {
        // 入力が2のべき乗であるかどうかを判断する関数
        
        return (n != 0) && ((n & (n - 1)) == 0)
    }
    
    func fft(_ x: [Float]) -> [Complex] {
        // 入力データの長さが2のべき乗であることを確認
        let log2n = vDSP_Length(log2(Double(x.count)).rounded(.up))
        let n = 1 << log2n
        
        // 入力データを複素数形式に変換
        var real = x
        var imaginary = [Float](repeating: 0.0, count: n)
        var splitComplexInput = DSPSplitComplex(realp: &real, imagp: &imaginary)
        
        // FFTの設定を行います。
        let fftSetup = vDSP_create_fftsetup(log2n, FFTRadix(kFFTRadix2))!
        
        // FFTを適用します。
        vDSP_fft_zip(fftSetup, &splitComplexInput, 1, log2n, FFTDirection(kFFTDirection_Forward))
        
        // 複素数の出力結果を作成します。
        var output = [Complex]()
        for i in 0..<n {
            output.append(Complex(re: splitComplexInput.realp[i], im: splitComplexInput.imagp[i]))
        }
        
        // FFT設定を解放します。
        vDSP_destroy_fftsetup(fftSetup)
        
        return output
    }
    
    var body: some View {
        // UI
        
        VStack{
            Chart {
                // データ構造からx, y値を取得してプロット
                
                ForEach(data) { shape in
                    LineMark(
                        x: .value("x", shape.xValue),
                        y: .value("y", shape.yValue)
                    )
                }
            }
            .padding(.all, 10)
            
            Button("Button"){
                // ボタンイベント
                
                // 関数を実行してx, y演算結果を取得
                let result = wavePoints(N: 8192, length: 10)
                let x = result.xArray
                let y = result.yArray
                // データ点数が2のべき乗であることを確認
                guard isPowerOfTwo(x.count) else {
                    // データ点数が2のべき乗でない場合、アラートを表示
                    showingAlert = true
                    return
                }
                
                // FFTを計算
                let fftResult = fft(y)
                
                // FFTの結果から振幅を計算
                var amp: [Float] = []
                for i in 0..<fftResult.count {
                    let re = fftResult[i].re
                    let im = fftResult[i].im
                    let magnitude = sqrt(re * re + im * im)
                    if i == 0 {
                        // 直流成分
                        amp.append(magnitude / Float(fftResult.count))
                    } else {
                        // 変動成分
                        amp.append(2 * magnitude / Float(fftResult.count))
                    }
                }
                
                // 周波数軸を計算
                let dt = x[1] - x[0] // 時間刻み
                var freq: [Float] = []
                for k in 0..<x.count {
                    freq.append(Float(k) / (Float(x.count) * dt))
                }
                
                // 周波数の上限を指定
                let freqLimit: Float = 300.0
                let limitIndex = Int(freqLimit * Float(fftResult.count) * dt)
                
                let filteredFreq = freq[..<limitIndex]
                let filteredAmp = amp[..<limitIndex]
                
                // プロットデータの全削除
                data.removeAll()
                
                // プロットデータの追加
                for i in 0..<filteredFreq.count {
                    data.append(.init(xValue: filteredFreq[i], yValue: filteredAmp[i]))
                }
                
            }
            .alert(isPresented: $showingAlert) {
                            Alert(title: Text("Error"), message: Text("データ点数が2の冪乗数ではありません。"), dismissButton: .default(Text("了解")))
                        }
        }
    }
}

struct ContentView_Previews: PreviewProvider {
    static var previews: some View {
        // プレビュー
        ContentView()
    }
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

import SwiftUI

import Charts

import Accelerate

struct PointsData: Identifiable {

// 点群データの構造体

var xValue: Float

var yValue: Float

var id = UUID()

}

struct ContentView: View {

// データを定義

@State private var data: [PointsData] = []

@State private var x: [Float] = []

@State private var y: [Float] = []

@State private var showingAlert = false

func linspace(start: Float, stop: Float, num: Int) -> [Float] {

// データ点数を指定して配列を作成する関数

let increment = (stop - start) / Float(num - 1)

return (0..<num).map { Float($0) * increment + start }

}

func wavePoints(N: Int, length: Float) -> (xArray: [Float], yArray: [Float]){

// プロットするデータを演算する関数

// データ点数Nと時間長lengthからx軸を作成

let xArray = linspace(start: 0.0, stop: length, num: N)

// 正弦波を作成。A:振幅、f:周波数[Hz]、dc:直流成分

let A: Float = 1

let f: Float = 100

let dc: Float = 1

let yArray = xArray.map{ A * sin(2 * .pi * f * $0) + dc }

return (xArray, yArray)

}

struct Complex {

// 複素数型の構造体を定義

// re:実部、im:虚部

var re: Float

var im: Float

// Complex + Complexの場合の演算処理

static func +(lhs: Complex, rhs: Complex) -> Complex {

return Complex(re: lhs.re + rhs.re, im: lhs.im + rhs.im)

}

// Complex * Complexの場合の演算処理

static func *(lhs: Complex, rhs: Complex) -> Complex {

return Complex(re: lhs.re * rhs.re - lhs.im * rhs.im, im: lhs.re * rhs.im + lhs.im * rhs.re)

}

// Complex * Floatの場合の演算処理

static func *(lhs: Complex, rhs: Float) -> Complex {

return Complex(re: lhs.re * rhs, im: lhs.im * rhs)

}

func isPowerOfTwo(_ n: Int) -> Bool {

// 入力が2のべき乗であるかどうかを判断する関数

return (n != 0) && ((n & (n - 1)) == 0)

}

func fft(_ x: [Float]) -> [Complex] {

// 入力データの長さが2のべき乗であることを確認

let log2n = vDSP_Length(log2(Double(x.count)).rounded(.up))

let n = 1 << log2n

// 入力データを複素数形式に変換

var real = x

var imaginary = [Float](repeating: 0.0, count: n)

var splitComplexInput = DSPSplitComplex(realp: &real, imagp: &imaginary)

// FFTの設定を行います。

let fftSetup = vDSP_create_fftsetup(log2n, FFTRadix(kFFTRadix2))!

// FFTを適用します。

vDSP_fft_zip(fftSetup, &splitComplexInput, 1, log2n, FFTDirection(kFFTDirection_Forward))

// 複素数の出力結果を作成します。

var output = [Complex]()

for i in 0..<n {

output.append(Complex(re: splitComplexInput.realp[i], im: splitComplexInput.imagp[i]))

}

// FFT設定を解放します。

vDSP_destroy_fftsetup(fftSetup)

return output

}

var body: some View {

// UI

VStack{

Chart {

// データ構造からx, y値を取得してプロット

ForEach(data) { shape in

LineMark(

x: .value("x", shape.xValue),

y: .value("y", shape.yValue)

)

}

.padding(.all, 10)

Button("Button"){

// ボタンイベント

// 関数を実行してx, y演算結果を取得

let result = wavePoints(N: 8192, length: 10)

let x = result.xArray

let y = result.yArray

// データ点数が2のべき乗であることを確認

guard isPowerOfTwo(x.count) else {

// データ点数が2のべき乗でない場合、アラートを表示

showingAlert = true

return

}

// FFTを計算

let fftResult = fft(y)

// FFTの結果から振幅を計算

var amp: [Float] = []

for i in 0..<fftResult.count {

let re = fftResult[i].re

let im = fftResult[i].im

let magnitude = sqrt(re * re + im * im)

if i == 0 {

// 直流成分

amp.append(magnitude / Float(fftResult.count))

} else {

// 変動成分

amp.append(2 * magnitude / Float(fftResult.count))

}

// 周波数軸を計算

let dt = x[1] - x[0] // 時間刻み

var freq: [Float] = []

for k in 0..<x.count {

freq.append(Float(k) / (Float(x.count) * dt))

}

// 周波数の上限を指定

let freqLimit: Float = 300.0

let limitIndex = Int(freqLimit * Float(fftResult.count) * dt)

let filteredFreq = freq[..<limitIndex]

let filteredAmp = amp[..<limitIndex]

// プロットデータの全削除

data.removeAll()

// プロットデータの追加

for i in 0..<filteredFreq.count {

data.append(.init(xValue: filteredFreq[i], yValue: filteredAmp[i]))

}

.alert(isPresented: $showingAlert) {

Alert(title: Text("Error"), message: Text("データ点数が2の冪乗数ではありません。"), dismissButton: .default(Text("了解")))

}

struct ContentView_Previews: PreviewProvider {

static var previews: some View {

// プレビュー

ContentView()

}

先ほどもデータ数が2の冪乗数でない時の処理はしていませんが、vDSPは自動的にゼロパディングによる自由なデータ数の処理が自動的に有効になっています。しかし振幅計算時に振幅補正係数を自分で書かないといけません。今回はSwiftになれていない筆者は少し面倒だったので、2の冪乗数でないデータが渡された場合に.alert()で警告文を出し計算しないという手段をとりました。

実行結果

こちらが実行結果です。

データ点数（let result = wavePoints(N: 8191, length: 10)）とNを2の冪乗数でない数を指定すると、以下の図に示すメッセージを出しています。

まとめ

今回はここまで！
Swiftを使ってフーリエ変換するコードを、まずはCooley-Tukey法をフルスクラッチすることで実装してみました。一度Pythonで検証していたといっても、やはりSwiftだと慣れていないので時間がかかりますね。

そして次にAccelerateというライブラリをimportし、vDSPメソッドでFFTを使ってみました。ライブラリを使うことで少し速度UPしたように感じました。
次回は録音再生アプリにFFTを実装する部分を書いてみたいと思います。

SwiftでもFFTをかけることができるようになりました！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

Swiftで高速フーリエ変換（FFT）を実装する方法

この記事の目標

Cooley-Tukey法によるFFTコードを自分で書いてみる

AccelerateでFFTを使ってみる

動作環境

SwiftでCooley-Tukey法によるFFTをフルスクラッチするコード

事前準備

データ点数を指定して配列を作成する関数

サンプル波形を生成する関数

複素数構造体を定義

FFT関数（Cooley-Tukey法）

全コード（コピペ動作OK）

実行結果

SwiftのAccelerate/vDSPでFFTをするコード

全コード

実行結果

まとめ

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

この記事の目標

Cooley-Tukey法によるFFTコードを自分で書いてみる

AccelerateでFFTを使ってみる

動作環境

SwiftでCooley-Tukey法によるFFTをフルスクラッチするコード

事前準備

データ点数を指定して配列を作成する関数

サンプル波形を生成する関数

複素数構造体を定義

FFT関数（Cooley-Tukey法）

全コード（コピペ動作OK）

実行結果

SwiftのAccelerate/vDSPでFFTをするコード

全コード

実行結果

まとめ

SNSでもご購読できます。

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

コメントを残すコメントをキャンセル