Strange results from Accelerate DFT

Question

Created 2w

Replies 4

Boosts 0

Participants 2

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

Answered by hilltopfarmer in 817045022

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
}

Boost

Answer 1

endecotp OP

2w

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.

According to Wikipedia:

https://en.wikipedia.org/wiki/Discrete_cosine_transform#Formal_definition

”This transform is exactly equivalent .. to a DFT of 4N real inputs of even symmetry where the even-indexed elements are zero. … DCT-II transformation is also possible using 2N signal followed by a multiplication by half shift”

What is the “multiplication by half shift” that they refer to? Is that what you’re missing? Have you tried the version with 4N inputs to the DFT?

1

Answer 2

hilltopfarmer OP

2w

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

0

Answer 3

hilltopfarmer OP

2w

Accepted Answer

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
}

0

Answer 4

hilltopfarmer OP

2w

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

0