Swift 通过线性调频 Z 变换 (CZT) 的反向 FFT (IFFT)
Swift Inverse FFT (IFFT) Via Chirp Z-Transform (CZT)
对于任意样本大小（样本不等于 2^N），我已经能够使用 iOS Accelerate 的 FFT 函数（仅适用于样本等于 2^N)。
结果很好，与任意长度序列（信号）的 Matlab FFT 输出相匹配。我粘贴下面的代码。
我的下一个挑战是使用 iOS Accelerate 的 FFT 函数（仅适用于等于 2^N 的样本）来完成任意样本大小（样本不等于 2^N）的逆 FFT。
由于我的 CZT 现在可以完成任意长度的 FFT（见下文），我希望逆 CZT (ICZT) 可以使用 iOS Accelerate 的 FFT 函数（仅适用于样本相等的样本）完成任意长度的 IFFT到 2^N).
任何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({ [=11=] / Double(x.count) }) // Normalize IFFT by sample count
imaginary = imaginary.map({ [=11=] / 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({ [=11=] + Double(m) })
var NM: [Double] = Array(stride(from: Double(-(m - 1)), through: Double(m - 1), by: 1.0))
NM = NM.map({ [=11=] + Double(m) })
var M: [Double] = Array(stride(from: Double(0), through: Double(m - 1), by: 1.0))
M = M.map({ [=11=] + 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({ [=11=] / 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: [=11=], imaginary: ) } // complexN 1
let fw = zip(fwr, fwi).map { Complex<Double>(real: [=11=], imaginary: ) } // complexN 2
let cc = zip(fg, fw).map { [=11=] * } // multiply above 2 complex numbers fg * fw
var ccR = cc.map { [=11=].real } // real part (vector) of complex multiply
var ccI = cc.map { [=11=].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: [=11=], imaginary: ) }
let www = zip(W2r, W2i).map { Complex<Double>(real: [=11=], imaginary: ) }
let y = zip(ggg, www).map { [=11=] * }
let yR = y.map { [=11=].real } // FFT real part (output vector)
let yI = y.map { [=11=].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
发布我的评论作为结束此问题的答案—
如果您确定要使用 ICZT 作为 IFFT 的等价物，那么让您的 dft 函数像您的 fft 一样接受 type: String 参数。当type为ifft时，你只需要翻转这里的符号：
WI.append(sin(-p1[i] * 2.0 * M_PI / Double(m)))
正向 FFT 为负，反向 FFT (IFFT) 为正。
这是我为演示 CZT 编写的一些 Octave/Matlab 代码：gist.github.com/fasiha/42a21405de92ea46f59e。该演示展示了如何使用 czt2 来执行 fft。 czt2 的第三个参数（在代码中称为 w）是 FFT 的 exp(-2j * pi / Nup)。只需将它与 exp(+2j * pi / Nup) 共轭即可获得 IFFT。
这就是翻转 WI 中 sin 中的符号的作用。
对于任意样本大小（样本不等于 2^N），我已经能够使用 iOS Accelerate 的 FFT 函数（仅适用于样本等于 2^N)。
结果很好，与任意长度序列（信号）的 Matlab FFT 输出相匹配。我粘贴下面的代码。
我的下一个挑战是使用 iOS Accelerate 的 FFT 函数（仅适用于等于 2^N 的样本）来完成任意样本大小（样本不等于 2^N）的逆 FFT。
由于我的 CZT 现在可以完成任意长度的 FFT（见下文），我希望逆 CZT (ICZT) 可以使用 iOS Accelerate 的 FFT 函数（仅适用于样本相等的样本）完成任意长度的 IFFT到 2^N).
任何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({ [=11=] / Double(x.count) }) // Normalize IFFT by sample count
imaginary = imaginary.map({ [=11=] / 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({ [=11=] + Double(m) })
var NM: [Double] = Array(stride(from: Double(-(m - 1)), through: Double(m - 1), by: 1.0))
NM = NM.map({ [=11=] + Double(m) })
var M: [Double] = Array(stride(from: Double(0), through: Double(m - 1), by: 1.0))
M = M.map({ [=11=] + 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({ [=11=] / 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: [=11=], imaginary: ) } // complexN 1
let fw = zip(fwr, fwi).map { Complex<Double>(real: [=11=], imaginary: ) } // complexN 2
let cc = zip(fg, fw).map { [=11=] * } // multiply above 2 complex numbers fg * fw
var ccR = cc.map { [=11=].real } // real part (vector) of complex multiply
var ccI = cc.map { [=11=].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: [=11=], imaginary: ) }
let www = zip(W2r, W2i).map { Complex<Double>(real: [=11=], imaginary: ) }
let y = zip(ggg, www).map { [=11=] * }
let yR = y.map { [=11=].real } // FFT real part (output vector)
let yI = y.map { [=11=].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
发布我的评论作为结束此问题的答案—
如果您确定要使用 ICZT 作为 IFFT 的等价物，那么让您的 dft 函数像您的 fft 一样接受 type: String 参数。当type为ifft时，你只需要翻转这里的符号：
WI.append(sin(-p1[i] * 2.0 * M_PI / Double(m)))
正向 FFT 为负，反向 FFT (IFFT) 为正。
这是我为演示 CZT 编写的一些 Octave/Matlab 代码：gist.github.com/fasiha/42a21405de92ea46f59e。该演示展示了如何使用 czt2 来执行 fft。 czt2 的第三个参数（在代码中称为 w）是 FFT 的 exp(-2j * pi / Nup)。只需将它与 exp(+2j * pi / Nup) 共轭即可获得 IFFT。
这就是翻转 WI 中 sin 中的符号的作用。