Swift Inverse FFT (IFFT) Via Chirp Z-Transfrom (CZT)

Question

Для произвольных размеров выборки (образцы, не равные 2^N), я смог реализовать БПФ через Z-преобразование чирпа (CZT) с помощью iOS Ускорить функцию FFT (которая работает только для образцов, равных 2^N).

Результаты хорошие и соответствуют выходному сигналу Matlab FFT для любой произвольной последовательности длин (сигнал). Я вставляю код ниже.

Следующая задача - использовать функцию FFT iOS Accelerate (которая работает только для образцов, равных 2^N) для выполнения обратного БПФ при произвольных размерах выборки (образцы не равны 2^N).

Поскольку мой CZT выполняет произвольную длину FFT теперь (см. Ниже), я надеюсь, что обратный CZT (ICZT) выполнит произвольную длину IFFT с использованием функции FFT iOS Accelerate (которая работает только для образцов, равных 2^N) ,

Любые предложения/рекомендации?

// 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 

Answer 1

проводки мой комментарий в ответ, чтобы закрыть этот вопрос-

Если вы уверены, что вы хотите использовать ICZT как эквивалент IFFT, а затем сделать ваш dft функции принимают type: String аргумент, как ваш fft. Когда typeifft это, все, что вам нужно, чтобы перевернуть знак здесь:

WI.append(sin(-p1[i] * 2.0 * M_PI/Double(m)))

Оставьте отрицательный для прямого преобразования Фурье, и положительное для обратного преобразования Фурье (БПФ).

---

Вот некоторые Octave/Matlab код, который я написал, чтобы продемонстрировать CZT: gist.github.com/fasiha/42a21405de92ea46f59e. Демонстрация показывает, как использовать czt2, чтобы сделать fft. Третий аргумент czt2 (в коде) называется w) - exp(-2j * pi/Nup) для БПФ. Просто соедините его с exp(+2j * pi/Nup), чтобы получить IFFT.

Это то, что переворачивает знак в sin в WI.