Skip to content

Instantly share code, notes, and snippets.

@moutend
Last active March 13, 2021 00:24
Show Gist options
  • Save moutend/065ed56bd45086d7a0f5c94916d602db to your computer and use it in GitHub Desktop.
Save moutend/065ed56bd45086d7a0f5c94916d602db to your computer and use it in GitHub Desktop.
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)
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment