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:
The dotted is the source, the solid is the output of IFFT(FFT(source)). Any ideas on where I'm going wrong?
