Post

Replies

Boosts

Views

Activity

Strange results from Accelerate DFT
Not sure if this is the right forum (don't see one for Maths...). There's no double-precision version of the newest Swift API for doing a DCT, so I tried to roll my own. Was expecting to have to scale the output, but by a constant factor... The results are out by a factor of about 0.2 at the start, and this falls smoothly to about 0.0067 by the end of the "period". I've included a direct calculation of the normalised DCT, as this was the kind of scaling error I was expecting to see. I think I'm doing it right: mirror the data into an array twice the size of the original, pass that in as the real part, pass zeroes in as the imaginary part, then take the real result. It's alternating zeroes all right, and it's in sync, but ... DO I need to scale the input? Chebyshev nodes. I'm at a loss. Anybody know what I'm doing wrong? PS Sorry for the huge long post. // Unnormalized values from SciPy, to save rewrites. // scipy.fft.dct([x for x in range(1,49)], type=2, norm=None) let fromSciPy = [ 2.35200000e+03, ... -3.27541474e-02] // Direct calculation which produces scaled output that agrees with SciPy // (although the zeroes don't seem to be exactly zero...) func dctII(input: [Double]) -> [Double] { let N = input.count let n = 1 / (2 * Double(N)) let factor = sqrt(2.0 / Double(N)) var result = (0..<N).map { k in factor * (0..<N).reduce(into: Double()) { sum_k, j in sum_k += input[j] * cos(Double.pi * Double(k) * (2 * Double(j) + 1) * n) } } result[0] /= sqrt(2.0) return result } let format = FloatingPointFormatStyle<Double>()... let factor = FloatingPointFormatStyle<Double>()... // According to the docs, an acceptable size. // .init() would fail if it wasn't. Effect's the same with monger sequences. let N = 48 let input = (1...N).map { Double($0) } // The bit that doesn't work. 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 to enforce even symmetry var real = [Double](repeating: 0, count: N * 2) var imag = [Double](repeating: 0, count: N * 2) 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) // Extract DCT-II coefficients from the real part of the first N terms var dctCoefficients = [Double](repeating: 0, count: N) for k in 0..<N { dctCoefficients[k] = re[k] } // Normalize to match SciPy's orthogonal normalization let scaleFactor = sqrt(2 / Double(N)) dctCoefficients[0] *= scaleFactor for k in 1..<N { dctCoefficients[k] *= scaleFactor } return dctCoefficients } let viaFFT = dctII_viaDFT(input: input) let direct = dctII(input: input) print(" SciPy Direct Drct/SciPy DFT DFT/SciPy FT/Direct") for i in direct.indices where viaFFT[i] != 0 { let truth = fromSciPy[i] let naive = direct[i] let weird = viaFFT[i] print("\(truth.formatted(format))\t\(naive.formatted(format))\t\((truth / naive).formatted(factor))\t\(weird.formatted(format))\t\t\((weird / truth).formatted(factor))\t\((weird / naive).formatted(factor))") } And here are the results (I've missed out the zero terms): SciPy Direct Drct/SciPy DFT DFT/SciPy FT/Direct +352.000000 +169.740979 3.856406 +480.099990 0.204124 2.828427 -933.609299 -095.286100 9.797959 -190.470165 0.204015 1.998929 -103.585662 -010.572167 9.797959 -021.042519 0.203141 1.990369 -037.182805 -003.794954 9.797959 -007.488532 0.201398 1.973287 -018.886898 -001.927636 9.797959 -003.754560 0.198792 1.947754 -011.356294 -001.159047 9.797959 -002.218278 0.195335 1.913881 -007.542756 -000.769829 9.797959 -001.440976 0.191041 1.871812 -005.347733 -000.545801 9.797959 -000.994300 0.185929 1.821728 -003.968777 -000.405062 9.797959 -000.714465 0.180021 1.763843 -003.045311 -000.310811 9.797959 -000.527882 0.173343 1.698404 -002.395797 -000.244520 9.797959 -000.397515 0.165922 1.625693 -001.920738 -000.196035 9.797959 -000.303073 0.157790 1.546021 -001.561880 -000.159409 9.797959 -000.232693 0.148983 1.459728 -001.283256 -000.130972 9.797959 -000.179063 0.139538 1.367185 -001.061666 -000.108356 9.797959 -000.137480 0.129495 1.268787 -000.881581 -000.089976 9.797959 -000.104818 0.118898 1.164955 -000.732264 -000.074736 9.797959 -000.078932 0.107791 1.056136 -000.606076 -000.061857 9.797959 -000.058319 0.096223 0.942793 -000.497433 -000.050769 9.797959 -000.041906 0.084243 0.825414 -000.402149 -000.041044 9.797959 -000.028916 0.071903 0.704500 -000.316996 -000.032353 9.797959 -000.018783 0.059254 0.580569 -000.239422 -000.024436 9.797959 -000.011098 0.046352 0.454153 -000.167336 -000.017079 9.797959 -000.005564 0.033251 0.325791 -000.098968 -000.010101 9.797959 -000.001980 0.020008 0.196034 -000.032754 -000.003343 9.797959 -000.000219 0.006679 0.065438
4
0
167
2w