Swift & v_DSP: Inverse FFT of FFT not returning approximately the same function

69 Views Asked by At

Doing some filtering in Swift for a signal processing application. While following the documentation on v_DSP, I put into practice carrying out an FFT() and then a IFFT() of the output. The implementation is as follows:

func filter() throws -> [Double] {
    let signal = self.map { Float($0) }
    let n = vDSP_Length(self.count)
    let log2n = vDSP_Length(log2(Float(n)))
    
    guard let fftSetUp = vDSP.FFT(
        log2n: log2n,
        radix: .radix2,
        ofType: DSPSplitComplex.self
    ) else {
        throw FilterError.fftUnavailable
    }
    
    let halfN = Int(n / 2)
    
    var forwardInputReal = [Float](repeating: 0,
                                   count: halfN)
    var forwardInputImag = [Float](repeating: 0,
                                   count: halfN)
    var forwardOutputReal = [Float](repeating: 0,
                                    count: halfN)
    var forwardOutputImag = [Float](repeating: 0,
                                    count: halfN)
    
    forwardInputReal.withUnsafeMutableBufferPointer { forwardInputRealPtr in
        forwardInputImag.withUnsafeMutableBufferPointer { forwardInputImagPtr in
            forwardOutputReal.withUnsafeMutableBufferPointer { forwardOutputRealPtr in
                forwardOutputImag.withUnsafeMutableBufferPointer { forwardOutputImagPtr in
                    
                    // Create a `DSPSplitComplex` to contain the signal.
                    var forwardInput = DSPSplitComplex(realp: forwardInputRealPtr.baseAddress!,
                                                       imagp: forwardInputImagPtr.baseAddress!)
                    
                    // Convert the real values in `signal` to complex numbers.
                    signal.withUnsafeBytes {
                        vDSP.convert(interleavedComplexVector: [DSPComplex]($0.bindMemory(to: DSPComplex.self)),
                                     toSplitComplexVector: &forwardInput)
                    }
                    
                    // Create a `DSPSplitComplex` to receive the FFT result.
                    var forwardOutput = DSPSplitComplex(realp: forwardOutputRealPtr.baseAddress!,
                                                        imagp: forwardOutputImagPtr.baseAddress!)
                    
                    // Perform the forward FFT.
                    fftSetUp.forward(input: forwardInput,
                                     output: &forwardOutput)
                }
            }
        }
    }
    
    var inverseOutputReal = [Float](repeating: 0,
                                    count: halfN)
    var inverseOutputImag = [Float](repeating: 0,
                                    count: halfN)


    let recreatedSignal: [Float] = forwardOutputReal.withUnsafeMutableBufferPointer { forwardOutputRealPtr in
        forwardOutputImag.withUnsafeMutableBufferPointer { forwardOutputImagPtr in
            inverseOutputReal.withUnsafeMutableBufferPointer { inverseOutputRealPtr in
                inverseOutputImag.withUnsafeMutableBufferPointer { inverseOutputImagPtr in
                    
                    // Create a `DSPSplitComplex` that contains the frequency-domain data.
                    let forwardOutput = DSPSplitComplex(realp: forwardOutputRealPtr.baseAddress!,
                                                        imagp: forwardOutputImagPtr.baseAddress!)
                    
                    // Create a `DSPSplitComplex` structure to receive the FFT result.
                    var inverseOutput = DSPSplitComplex(realp: inverseOutputRealPtr.baseAddress!,
                                                        imagp: inverseOutputImagPtr.baseAddress!)
                    
                    // Perform the inverse FFT.
                    fftSetUp.inverse(input: forwardOutput,
                                     output: &inverseOutput)
                    
                    // Return an array of real values from the FFT result.
                    let scale = 1 / Float(n * 2)
                    return [Float](fromSplitComplex: inverseOutput,
                                   scale: scale,
                                   count: Int(n))
                }
            }
        }
    }
    
    return recreatedSignal.map { Double($0) }
}

However when I plot the source and the output of the IFFT(FFT(source)), the data seems to be off by a scaling factor, and constant near the end:

Source vs. IFFT(FFT(source))

The dotted is the source, the solid is the output of IFFT(FFT(source)). Any ideas on where I'm going wrong?

0

There are 0 best solutions below