Swift vDSP rFFT not same as Python np.fft.rfft()

Question

I am trying to implement real FFT in iOS, for that i am using Accelerate Framework. Here is my Code for Swift.

class FFT {

private var fftSetup: FFTSetup?
private var log2n: Float?
private var length: Int?

func initialize(count: Int){
    length = count
    log2n = log2(Float(length!))
    self.fftSetup = vDSP_create_fftsetup(vDSP_Length(log2n!), FFTRadix(kFFTRadix2))!
}

func computeFFT(input: [Float]) -> ([Float], [Float]) {
    
    var real = input
    var imag = [Float](repeating: 0.0, count: input.count)
    var splitComplexBuffer = DSPSplitComplex(realp: &real, imagp: &imag)

    let halfLength = (input.count/2) + 1
    real = [Float](repeating: 0.0, count: halfLength)
    imag = [Float](repeating: 0.0, count: halfLength)

    // input is alternated across the real and imaginary arrays of the DSPSplitComplex structure
    splitComplexBuffer = DSPSplitComplex(fromInputArray: input, realParts: &real, imaginaryParts: &imag)

    // even though there are 2 real and 2 imaginary output elements, we still need to ask the fft to process 4 input samples
    vDSP_fft_zrip(fftSetup!, &splitComplexBuffer, 1, vDSP_Length(log2n!), FFTDirection(FFT_FORWARD))

    // zrip results are 2x the standard FFT and need to be scaled
    var scaleFactor = Float(1.0/2.0)
    vDSP_vsmul(splitComplexBuffer.realp, 1, &scaleFactor, splitComplexBuffer.realp, 1, vDSP_Length(halfLength))
    vDSP_vsmul(splitComplexBuffer.imagp, 1, &scaleFactor, splitComplexBuffer.imagp, 1, vDSP_Length(halfLength))
    
    return (real, imag)
    
}

func computeIFFT(real: [Float], imag: [Float]) -> [Float]{
    
    var real = [Float](real)
    var imag = [Float](imag)
    
    var result : [Float] = [Float](repeating: 0.0, count: length!)
    var resultAsComplex : UnsafeMutablePointer<DSPComplex>? = nil

    result.withUnsafeMutableBytes {
        resultAsComplex = $0.baseAddress?.bindMemory(to: DSPComplex.self, capacity: 512)
    }
    
    var splitComplexBuffer = DSPSplitComplex(realp: &real, imagp: &imag)
    
    vDSP_fft_zrip(fftSetup!, &splitComplexBuffer, 1, vDSP_Length(log2n!), FFTDirection(FFT_INVERSE));
    
    vDSP_ztoc(&splitComplexBuffer, 1, resultAsComplex!, 2, vDSP_Length(length! / 2));
    //
    //// Neither the forward nor inverse FFT does any scaling. Here we compensate for that.
    var scale : Float = 1.0/Float(length!);
    var copyOfResult = result;
    vDSP_vsmul(&result, 1, &scale, &copyOfResult, 1, vDSP_Length(length!));
    result = copyOfResult
    
    return result
}

func deinitialize(){
    vDSP_destroy_fftsetup(fftSetup)
}

}

Here is my code of Python for computing of rFFT and irFFT

# calculate fft of input block
    in_block_fft = np.fft.rfft(np.squeeze(in_buffer)).astype("complex64")

# apply mask and calculate the ifft
    estimated_block = np.fft.irfft(in_block_fft * out_mask)

Question:

Swift If i compute rFFT for a 512 frame and apply irFFT to the resultant of rFFT, I get the same original Array.

Python Same goes for python, if i take rFFT and irFFT i get the original Array in return.

Problem The problem occurs if i compare the results of Swift rFFT and Python rFFT. Their results are different in decimal values. Sometimes the real part is same but the imaginary is totally different.

I tried different Framework in Python like Numpy, SciPy and TensorFlow and their results are exactly the same (Slight different in decimal part). But when i calculate rfft on same input for iOS using the Swift code above, results are different.

If anybody who has experience with Accelerate Framework and holds the knowledge of FFT help me out with this one would be very helpful. I have limited knowledge of FFT.

Answer 1

Two observations: first your calls to DSPSplitComplex(realp: &real, imagp: &imag) create temporary pointers and the fromInputArray initializer is deprecated. Take a look at Data Packing for Fourier Transforms for best practice.

Secondly, the bindMemory capacity is the number of elements. The following lines create 512 complex elements where they should (if I understand correctly) create 256 complex elements to represent 512 real elements.

        result.withUnsafeMutableBytes {
            resultAsComplex = $0.baseAddress?.bindMemory(to: DSPComplex.self, capacity: 512)
        }

simon

Answer 2

Yes, it has differences.

I partially ported this functionality from python librosa to iOS with the same data.

https://github.com/dhrebeniuk/RosaKit