高速フーリエ変換を実施する(vDSPを使う)
GZYtechiPhone や Mac などでフーリエ変換を実施する場合は、Apple社の Accelerate フレームワークという信号処理用の高速なライブラリを使いましょう。そこに vDSP_fft_zip という FFT に使えるメソッドがあります。Apple開発者サイトを調べまくって、なんとか実行できるコードを作ることができましたので紹介します。
ちなみに下記の記事は、むかしむかしに独自で書いたC言語のFFTのコードを Swift に翻訳したものを紹介しています。( 独自で書いたというよりは、トラ技に乗っていたコードを自分テイストに書き換えただけですが... )
GZYtech 
自分でいろいろ高速化を考えたりするよりも、プライド捨ててサッサと vDSP を使ったほうが速いのがくやしいですッ！
import SwiftUI
// データ配列の要素をなす構造体の型.
struct X_FX {
var x: Double
var fx: Double
}
// スライダ用の定数.
struct SLIDER {
static let MIN: Double = 0.0
static let MAX: Double = 1022.0
static let STEP: Double = 1.0
}
struct ContentView: View {
// 背景の色と波形の色.
let COLOR_BG = Color( red: 0.0, green: 0.0, blue: 0.3, opacity: 1.0 )
let COLOR_WAVE = Color( red: 0.0, green: 1.0, blue: 0.0, opacity: 1.0 )
// データの配列.
let NUM_DATA = 2048
@State var Data: [X_FX] = []
// フーリエ変換オブジェクト.
let Fourier = MyFourierObject()
// スライダと関係するフィルタリングの次数.
@State var SliderValue: Double = 60.0
@State var Threshold = 60
// ローパスフィルタリングするかハイパスフィルタリングするか.
@State var ModePassLow: Bool = true
// f(x) 方向の描画倍率
let DRAW_RATE = 0.8
// FFTとIFFTの実行時間.
@State var MsecExec: Int = 0
var body: some View {
GeometryReader {
geom in
ZStack{
// 画面のうち ZStack が占めるサイズ.
let RectW = geom.size.width
let RectH = geom.size.height
// 波形描画の背景の矩形を描画する.
Rectangle()
.foregroundColor( COLOR_BG )
.frame( width: RectW, height: RectH )
// 目安になるセンター十字線と振幅の高さを描画する.
Path{ path in
let plot_x = RectW * 0.5
let plot_y = RectH * 0.5
// センターの垂直線.
let pnt00 = CGPoint( x: plot_x, y: 0 )
let pnt01 = CGPoint( x: plot_x, y: RectH )
path.move( to: pnt00 )
path.addLine( to: pnt01 )
// センターの水平線.
let pnt10 = CGPoint( x: 0, y: plot_y )
let pnt11 = CGPoint( x: RectW, y: plot_y )
path.move( to: pnt10 )
path.addLine( to: pnt11 )
// 振幅 f(x) のプラス側のめやす水平線.
let tmp0 = ( plot_y - ( plot_y * DRAW_RATE ))
let pnt20 = CGPoint( x: 0, y: tmp0 )
let pnt21 = CGPoint( x: RectW, y: tmp0 )
path.move( to: pnt20 )
path.addLine( to: pnt21 )
// 振幅 f(x) のプラス側のめやす水平線.
let tmp1 = ( plot_y + ( plot_y * DRAW_RATE ))
let pnt30 = CGPoint( x: 0, y: tmp1 )
let pnt31 = CGPoint( x: RectW, y: tmp1 )
path.move( to: pnt30 )
path.addLine( to: pnt31 )
}
.stroke( .gray, lineWidth: 1.0 )
.frame( width: RectW, height: RectH ) // ここが重要.
// 波形を描画する.
Path{ path in
let amp = RectH * 0.5 * DRAW_RATE
let off_y = RectH * 0.5
if ( Data.count > 0 )
{
var plot_x = Double(0)
var plot_y = -( Data[0].fx ) * amp + off_y
var pnt = CGPoint( x: plot_x, y: plot_y )
path.move( to: pnt )
let loop = Int( RectW )
for i in 1 ..< loop {
let rate = Double(i)/Double( RectW - 1 )
let n = Int( rate * Double( Data.count - 1 ))
plot_x = Double(i)
plot_y = -( Data[n].fx ) * amp + off_y
pnt = CGPoint( x: plot_x, y: plot_y )
path.addLine( to: pnt )
}
}
}
.stroke( COLOR_WAVE, lineWidth: 1.5 )
.frame( width: RectW, height: RectH ) // ここが重要.
VStack{
// データをセットするボタンを水平に並べる.
HStack{
// ローパスフィルタを選択するボタン.
Button( action: {
ModePassLow = true
processing( &Data, NUM_DATA, Threshold, ModePassLow )
}){
Text( "LowPass" )
.foregroundColor( .white )
}
.padding()
.background( .blue )
// ハイパスフィルタを選択するボタン.
Button( action: {
ModePassLow = false
processing( &Data, NUM_DATA, Threshold, ModePassLow )
}){
Text( "HighPass" )
.foregroundColor( .white )
}
.padding()
.background( .blue )
// フィルタ次数をひとつ増加するボタン.
Button( action: {
SliderValue += 1
if SliderValue > SLIDER.MAX
{
SliderValue = SLIDER.MAX
}
Threshold = Int( SliderValue )
processing( &Data, NUM_DATA, Threshold, ModePassLow )
}){
Text( "+" )
.foregroundColor( .white )
}
.padding()
.background( .blue )
// フィルタ次数をひとつ減らすボタン.
Button( action: {
SliderValue -= 1
if SliderValue < SLIDER.MIN
{
SliderValue = SLIDER.MIN
}
Threshold = Int( SliderValue )
processing( &Data, NUM_DATA, Threshold, ModePassLow )
}){
Text( "-" )
.foregroundColor( .white )
}
.padding()
.background( .blue )
}
// スライダの値を示す.
Text( String( format: "%f", SliderValue ) )
.foregroundColor( .white )
// フィルタ次数の値を示す.
Text( String( format: "%d / %d", Threshold, NUM_DATA ) )
.foregroundColor( .white )
// 実行時間を示す.
Text( String( format: "%d msec", MsecExec ) )
.foregroundColor( .white )
// フィルタ次数を決定するスライダ.
Slider(
value: $SliderValue,
in: SLIDER.MIN...SLIDER.MAX,
step: SLIDER.STEP,
onEditingChanged: {
editing in
// これをコメントインしてもなんだか違いがわからない.
//Threshold = Int( SliderValue )
//processing( &Data, NUM_DATA, Threshold, ModePassLow )
}
).onChange( of: SliderValue ){ new_value in
Threshold = Int( new_value )
processing( &Data, NUM_DATA, Threshold, ModePassLow )
}
}
.frame( maxWidth: .infinity, maxHeight: .infinity, alignment: .bottom )
.padding( .bottom )
}
.onAppear(){
// 起動時に矩形波データを仕込む.
SetDataSquare( &Data, NUM_DATA )
}
}
}
// フーリエ変換して周波数領域でフィルタリングして逆フーリエ変換する.
func processing( _ data: inout[X_FX], _ num_data: Int, _ threshold: Int, _ mode_filter: Bool )->Void
{
// 矩形波を仕込む.
let _ = SetDataSquare( &data, num_data )
// 実部と虚部のデータ領域を確保する.
let data_re = UnsafeMutablePointer<Float>.allocate( capacity: num_data )
let data_im = UnsafeMutablePointer<Float>.allocate( capacity: num_data )
let buf_for_dsp = UnsafeMutablePointer<Float>.allocate( capacity: num_data )
// f(x) のデータを実部に仕込み、虚部に 0.0 を仕込む.
for n in 0 ..< num_data
{
data_re[n] = Float( data[n].fx )
data_im[n] = Float( 0.0 )
}
// 時間計測を開始する.
let date_beg = Date()
// フーリエ変換で周波数領域の実部と虚部のデータを取得する.
Fourier.Transform(
MyFourierObject.TransformType.TimeToFreq,
data_re,
data_im,
buf_for_dsp,
num_data
)
// 周波数領域でフィルタリングを実施する.
if ( mode_filter )
{
// ローパスフィルタ.
let _ = FilteringPassLow(
data_re,
data_im,
num_data,
threshold
)
}
else
{
// ハイパスフィルタ.
let _ = FilteringPassHigh(
data_re,
data_im,
num_data,
threshold
)
}
// フーリエ逆変換する.
Fourier.Transform(
MyFourierObject.TransformType.FreqToTime,
data_re,
data_im,
buf_for_dsp,
num_data
)
// 時間計測を終わる.
let date_fin = Date()
// 実行時間を取得する.
let sec = date_fin.timeIntervalSince( date_beg )
let msec = sec * 1000.0
MsecExec = Int( msec.rounded())
// 実部のデータを元のデータ f(x) に書き戻す.
for n in 0 ..< num_data
{
data[n].fx = Double( data_re[n] )
}
data_re.deallocate()
data_im.deallocate()
buf_for_dsp.deallocate()
return
}
// サイン波形を格納する.
func SetDataSin( _ data: inout [X_FX], _ req_num_data: Int )->Void {
// いったん配列をクリアする.
data.removeAll()
let RANGE = 2.0
let MUL = 2.0
for n in 0 ..< req_num_data {
let rate = Double(n)/Double( req_num_data - 1)
var v = RANGE * rate
v -= 1.0 // この時点で -1.0 〜 0.0 〜 +1.0 になる.
let t = v * ( 2.0 * ( Double.pi )) * MUL
let ft = sin( t )
// データ要素を追加する.
let element = X_FX( x: t, fx: ft )
data.append( element )
}
}
// 矩形波を格納する.
func SetDataSquare( _ data: inout [X_FX], _ req_num_data: Int )->Void {
// いったん配列をクリアする.
data.removeAll()
// サイン波形をしこむ.
SetDataSin( &data, req_num_data )
// サイン波形の正負を評価して、矩形波をつくる.
for n in 0 ..< req_num_data {
if ( data[n].fx >= 0.0 )
{
data[n].fx = +1.0;
}
else
{
data[n].fx = -1.0
}
}
}
// ローパスフィルタリングを実施する.
func FilteringPassLow(
_ data_re: UnsafeMutablePointer<Float>,
_ data_im: UnsafeMutablePointer<Float>,
_ num_data: Int,
_ threshold: Int
)->Bool
{
let beg = threshold
let fin = num_data - threshold
if ( beg < 0 ){ return false }
if ( beg >= num_data ) { return false }
if ( fin < 0 ){ return false }
if ( fin >= num_data ) { return false }
for k in beg ... fin
{
data_re[k] = 0.0
data_im[k] = 0.0
}
return true
}
// ハイパスフィルタリングを実施する.
func FilteringPassHigh(
_ data_re: UnsafeMutablePointer<Float>,
_ data_im: UnsafeMutablePointer<Float>,
_ num_data: Int,
_ threshold: Int
)->Bool
{
let half = num_data >> 1
let beg = half - threshold
let fin = half + threshold
if ( beg < 0 ){ return false }
if ( beg >= num_data ) { return false }
if ( fin < 0 ){ return false }
if ( fin >= num_data ) { return false }
var a = 0
var b = 0
a = 0
b = beg
for k in a ..< b
{
data_re[k] = 0.0
data_im[k] = 0.0
}
a = fin + 1
b = num_data - 1
for k in a ..< b
{
data_re[k] = 0.0
data_im[k] = 0.0
}
return true
}
}
struct ContentView_Previews: PreviewProvider {
static var previews: some View {
ContentView()
}
}
下記が vDSP を用いたフーリエ変換とフーリエ逆変換を実行するクラスです。
Transform というメソッドの buffer_for_dsp というメモリ領域は、内部で DSP 乗算をするので、それに用いる領域です。このコードを改造して引数を減らしたい場合は、メソッドの内部で buffer_for_dsp メモリをアロケートしてください。中でいちいちメモリ確保するので、その分が速度に影響するかもしれませんが。
データの個数は 2 のべき乗で実行してください。
import Accelerate
class MyFourierObject {
public enum TransformType {
case TimeToFreq
case FreqToTime
}
init()
{
}
public func Transform(
_ transform_type: TransformType,
_ data_re: UnsafeMutablePointer<Float>,
_ data_im: UnsafeMutablePointer<Float>,
_ buffer_for_dsp: UnsafeMutablePointer<Float>,
_ num_data: Int
) -> Void
{
// 複素型データの準備をする.
var split_complex = DSPSplitComplex( realp: data_re, imagp: data_im )
// FFTの準備をする.
let dsp_length = vDSP_Length( log2( Float( num_data )))
let fft_setup = vDSP_create_fftsetup( dsp_length, FFTRadix( kFFTRadix2 ))
let MY_FFT_STRIDE = 1
// 場合分け.
if ( transform_type == TransformType.TimeToFreq )
{
// TimeToFreq.
vDSP_fft_zip(
fft_setup!,
&split_complex,
MY_FFT_STRIDE,
dsp_length,
FFTDirection( FFT_FORWARD )
)
}
else
{
// FreqToTime.
vDSP_fft_zip(
fft_setup!,
&split_complex,
MY_FFT_STRIDE,
dsp_length,
FFTDirection( FFT_INVERSE )
)
}
// FFTの後処理をする.
vDSP_destroy_fftsetup( fft_setup )
// 場合分け.
if ( transform_type == TransformType.TimeToFreq )
{
// NOP.
}
else
{
// スケーリングのための係数を計算する.
var scale = Float( 1.0 )/Float( num_data )
// DSP乗算のストライド.
let MY_DSP_STRIDE = 1
// 実部のDSP乗算を実施する.
vDSP_vsmul(
data_re,
MY_DSP_STRIDE,
&scale,
buffer_for_dsp,
MY_DSP_STRIDE,
vDSP_Length( num_data )
)
// 引数の実部データ列に戻す.
for n in 0 ..< num_data
{
data_re[n] = buffer_for_dsp[n]
}
// 虚部のDSP乗算を実施する.
vDSP_vsmul(
data_im,
MY_DSP_STRIDE,
&scale,
buffer_for_dsp,
MY_DSP_STRIDE,
vDSP_Length( num_data )
)
// 引数の虚部データ列に戻す.
for n in 0 ..< num_data
{
data_im[n] = buffer_for_dsp[n]
}
}
return
}
}
