Edit: Don't use the reduce formulation (typed in haste...) You've still got to skip the zeroes...

Sorted. Turns out you mustn't ignore the imaginary part... Here's what to do, in case anyone else has a "moment": func dctII_viaDFT(input: [Double]) -> [Double] { let N = input.count // Initialize DFT object for 2N points guard let dft = try? vDSP.DiscreteFourierTransform( count: N * 2, direction: .forward, transformType: .complexComplex, ofType: Double.self ) else { fatalError("Failed to create DFT object") } // Extend the input signal var real = [Double](repeating: 0, count: N * 2) // Imaginary part is all zeroes let imag = [Double](repeating: 0, count: N * 2) // Real part is even symmetric for i in 0..<N { real[i] = input[i] real[(N * 2) - 1 - i] = input[i] } // Compute the DFT let (re, im) = dft.transform(real: real, imaginary: imag) // Sort out the mess (including optional normalisation) let halfShift = .pi / Double(N * 2) // normalisation factor, which doesn't have to be here... let scaleFactor = 1 / sqrt(2 * Double(N)) var sine = 0.0 var cosine = 0.0 var result = re.indices.reduce(into: [Double]()) { coeffs, k in __sincos(halfShift * Double(k), &sine, &cosine) // Doh! // scaleFactor is only for normalised output let coeff = (re[k] * cosine + im[k] * sine) * scaleFactor coeffs.append(coeff) } // also only for normalised output result[0] /= sqrt(2) return result }

Thanks for the suggestions. The half-shift, in this case, would be cos((pi * k)/n) for k over the whole series. But the result doesn't look shifted, it just looks wrong. There may be something there, though. I'll definitely try interleaving the real and imaginary parts, though, because, as you say, that isn't supposed to need the half-shift. Thanks again