Pour les tailles d'échantillons arbitraires (échantillons non égal à 2^N), j'ai été en mesure de mettre en œuvre la FFT via le gazouillement transformée en Z (CZT) en utilisant iOS Fonction FFT d'Accelerate (qui ne fonctionne que pour des échantillons égaux à 2^N).Swift FFT inverse (IFFT) Via Chirp Z-Transfrom (CZT)
Les résultats sont bons et correspondent à la sortie Matlab FFT pour chaque séquence de longueur arbitraire (signal). Je colle le code ci-dessous. Mon prochain défi est d'utiliser la fonction FFT d'iOS Accelerate (qui ne fonctionne que pour des échantillons égaux à 2^N) pour réaliser une FFT inverse sur des tailles d'échantillons arbitraires (échantillons non égaux à 2^N). Puisque mon CZT accomplit maintenant la longueur arbitraire FFT (voir ci-dessous), j'espère qu'un CZT inverse (ICZT) accomplirait une longueur arbitraire IFFT en utilisant la fonction FFT d'iOS Accelerate (cela fonctionne seulement pour des échantillons égaux à 2^N) .
Toutes les suggestions/guidence?
// FFT IOS ACCELERATE FRAMEWORK (works only for 2^N samples)
import Accelerate
public func fft(x: [Double], y: [Double], type: String) -> ([Double], [Double]) {
var real = [Double](x)
var imaginary = [Double](y)
var splitComplex = DSPDoubleSplitComplex(realp: &real, imagp: &imaginary)
let length = vDSP_Length(floor(log2(Float(real.count))))
let radix = FFTRadix(kFFTRadix2)
let weights = vDSP_create_fftsetupD(length, radix)
switch type.lowercased() {
case ("fft"): // CASE FFT
vDSP_fft_zipD(weights!, &splitComplex, 1, length, FFTDirection(FFT_FORWARD))
vDSP_destroy_fftsetup(weights)
case ("ifft"): // CASE INVERSE FFT
vDSP_fft_zipD(weights!, &splitComplex, 1, length, FFTDirection(FFT_INVERSE))
vDSP_destroy_fftsetup(weights)
real = real.map({ $0/Double(x.count) }) // Normalize IFFT by sample count
imaginary = imaginary.map({ $0/Double(x.count) }) // Normalize IFFT by sample count
default: // DEFAULT CASE (FFT)
vDSP_fft_zipD(weights!, &splitComplex, 1, length, FFTDirection(FFT_FORWARD))
vDSP_destroy_fftsetup(weights)
}
return (real, imaginary)
}
// END FFT IOS ACCELERATE FRAMEWORK (works only for 2^N samples)
// DEFINE COMPLEX NUMBERS
struct Complex<T: FloatingPoint> {
let real: T
let imaginary: T
static func +(lhs: Complex<T>, rhs: Complex<T>) -> Complex<T> {
return Complex(real: lhs.real + rhs.real, imaginary: lhs.imaginary + rhs.imaginary)
}
static func -(lhs: Complex<T>, rhs: Complex<T>) -> Complex<T> {
return Complex(real: lhs.real - rhs.real, imaginary: lhs.imaginary - rhs.imaginary)
}
static func *(lhs: Complex<T>, rhs: Complex<T>) -> Complex<T> {
return Complex(real: lhs.real * rhs.real - lhs.imaginary * rhs.imaginary,
imaginary: lhs.imaginary * rhs.real + lhs.real * rhs.imaginary)
}
}
extension Complex: CustomStringConvertible {
var description: String {
switch (real, imaginary) {
case (_, 0):
return "\(real)"
case (0, _):
return "\(imaginary)i"
case (_, let b) where b < 0:
return "\(real) - \(abs(imaginary))i"
default:
return "\(real) + \(imaginary)i"
}
}
}
// DEFINE COMPLEX NUMBERS
// DFT BASED ON CHIRP Z TRANSFORM (CZT)
public func dft(x: [Double]) -> ([Double], [Double]) {
let m = x.count // number of samples
var N: [Double] = Array(stride(from: Double(0), through: Double(m - 1), by: 1.0))
N = N.map({ $0 + Double(m) })
var NM: [Double] = Array(stride(from: Double(-(m - 1)), through: Double(m - 1), by: 1.0))
NM = NM.map({ $0 + Double(m) })
var M: [Double] = Array(stride(from: Double(0), through: Double(m - 1), by: 1.0))
M = M.map({ $0 + Double(m) })
let nfft = Int(pow(2, ceil(log2(Double(m + m - 1))))) // fft pad
var p1: [Double] = Array(stride(from: Double(-(m - 1)), through: Double(m - 1), by: 1.0))
p1 = (zip(p1, p1).map(*)).map({ $0/Double(2) }) // W = WR + j*WI has to be raised to power p1
var WR = [Double]()
var WI = [Double]()
for i in 0 ..< p1.count { // Use De Moivre's formula to raise to power p1
WR.append(cos(p1[i] * 2.0 * M_PI/Double(m)))
WI.append(sin(-p1[i] * 2.0 * M_PI/Double(m)))
}
var aaR = [Double]()
var aaI = [Double]()
for j in 0 ..< N.count {
aaR.append(WR[Int(N[j] - 1)] * x[j])
aaI.append(WI[Int(N[j] - 1)] * x[j])
}
let la = nfft - aaR.count
let pad: [Double] = Array(repeating: 0, count: la) // 1st zero padding
aaR += pad
aaI += pad
let (fgr, fgi) = fft(x: aaR, y: aaI, type: "fft") // 1st FFT
var bbR = [Double]()
var bbI = [Double]()
for k in 0 ..< NM.count {
bbR.append((WR[Int(NM[k] - 1)])/(((WR[Int(NM[k] - 1)])) * ((WR[Int(NM[k] - 1)])) + ((WI[Int(NM[k] - 1)])) * ((WI[Int(NM[k] - 1)])))) // take reciprocal
bbI.append(-(WI[Int(NM[k] - 1)])/(((WR[Int(NM[k] - 1)])) * ((WR[Int(NM[k] - 1)])) + ((WI[Int(NM[k] - 1)])) * ((WI[Int(NM[k] - 1)])))) // take reciprocal
}
let lb = nfft - bbR.count
let pad2: [Double] = Array(repeating: 0, count: lb) // 2nd zero padding
bbR += pad2
bbI += pad2
let (fwr, fwi) = fft(x: bbR, y: bbI, type: "fft") // 2nd FFT
let fg = zip(fgr, fgi).map { Complex<Double>(real: $0, imaginary: $1) } // complexN 1
let fw = zip(fwr, fwi).map { Complex<Double>(real: $0, imaginary: $1) } // complexN 2
let cc = zip(fg, fw).map { $0 * $1 } // multiply above 2 complex numbers fg * fw
var ccR = cc.map { $0.real } // real part (vector) of complex multiply
var ccI = cc.map { $0.imaginary } // imag part (vector) of complex multiply
let lc = nfft - ccR.count
let pad3: [Double] = Array(repeating: 0, count: lc) // 3rd zero padding
ccR += pad3
ccI += pad3
let (ggr, ggi) = fft(x: ccR, y: ccI, type: "ifft") // 3rd FFT (IFFT)
var GGr = [Double]()
var GGi = [Double]()
var W2r = [Double]()
var W2i = [Double]()
for v in 0 ..< M.count {
GGr.append(ggr[Int(M[v] - 1)])
GGi.append(ggi[Int(M[v] - 1)])
W2r.append(WR[Int(M[v] - 1)])
W2i.append(WI[Int(M[v] - 1)])
}
let ggg = zip(GGr, GGi).map { Complex<Double>(real: $0, imaginary: $1) }
let www = zip(W2r, W2i).map { Complex<Double>(real: $0, imaginary: $1) }
let y = zip(ggg, www).map { $0 * $1 }
let yR = y.map { $0.real } // FFT real part (output vector)
let yI = y.map { $0.imaginary } // FFT imag part (output vector)
return (yR, yI)
}
// END DFT BASED ON CHIRP Z TRANSFORM (CZT)
// CHIRP DFT (CZT) TEST
let x: [Double] = [1, 2, 3, 4, 5] // arbitrary sample size
let (fftR, fftI) = dft(x: x)
print("DFT Real Part:", fftR)
print(" ")
print("DFT Imag Part:", fftI)
// Matches Matlab FFT Output
// DFT Real Part: [15.0, -2.5000000000000018, -2.5000000000000013, -2.4999999999999991, -2.499999999999996]
// DFT Imag Part: [-1.1102230246251565e-16, 3.4409548011779334, 0.81229924058226477, -0.81229924058226599, -3.4409548011779356]
// END CHIRP DFT (CZT) TEST
Oh mon Dieu, tant de questions. Utilisez-vous simplement CZT pour obtenir des FFT non puissantes sur deux en utilisant la FFT d'Accelerate (qui ne fonctionne que pour la puissance de deux)? Pourquoi n'utilisez-vous pas une autre bibliothèque FFT non contrainte, comme FFTW ou une autre bibliothèque Apple? Toute accélération obtenue en utilisant Accelerate sous le capot sera probablement complètement évitée par l'énorme charge de calcul de la CZT (mise en place de toutes ces exponentielles). –
Si vous voulez vraiment vraiment vraiment vraiment utiliser ICZT pour obtenir des IFFT non-power-of-two de cette façon, faites que votre fonction 'dft' accepte un argument' type: String' comme votre 'fft'. Quand 'type' est' ifft', je pense que tout ce dont vous avez besoin est de retourner le signe ici: 'WI.append (sin (-p1 [i] * 2.0 * M_PI/Double (m)))', ont un négatif pour l'avant , et positif pour inverse. Je * pense * que c'est tout ce que vous devez faire pour changer la direction de CZT → ICZT. –
Voir ce code Matlab/Octave: https://gist.github.com/fasiha/42a21405de92ea46f59e La démo montre comment utiliser '' czt2' faire fft' avec le troisième argument à 'czt2', appelé' W', étant 'exp (-2j * pi/Nup)'. Pour le faire correspondre 'ifft', changez simplement le troisième argument en' exp (+ 2j * pi/Nup) ', c'est-à-dire, conjuguez' w'. C'est ce que fait le signe dans l'exposant pour 'WI '. –