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