1

I try to transform a simple cosine signal using the Discrete Cosine Transform (DCT) scipy.fft.dct, however it seems there is an issue as there is power in frequencies that should not exist. Suppose a domain from zero to one, both endpoints included, for the cosine function:

import numpy as np
x = np.linspace(0, 1, 8, endpoint = True)
f = np.cos(1 * np.pi * x)

This simple signal offers a single frequency, so I do expect significant powers only at a single frequency of the DCT:

import scipy.fft
f_FT = scipy.fft.dct(f, type = 1, norm = "ortho")

I select the DCT type I according to the Wikipedia classification (that is also referenced in SciPy's documentation) because the endpoints are included and the signal is even at both boundaries. But this yields as result:

array([ 3.35699888e-16,  2.09223516e+00, -1.48359792e-17,  2.21406462e-01,
       -1.92867730e-16,  2.21406462e-01,  1.18687834e-16,  1.56558011e-01])

Thus, there is still significant energy in k=3pi, 5pi, 7pi (second and last column).

Am I doing something wrong? As written above, I expect only power at k=1pi. The Discrete Sine Transform (DST) does not offer this kind of problem - there, I find only power in frequencies that I generate.

Thank you in advance for your help.

janw
  • 8,758
  • 11
  • 40
  • 62
phil
  • 11
  • 2

1 Answers1

0

Update - origin of problem found

I found that the problem is caused by norm = "ortho". Apparently, the library modifies the first and last point of the input signal before the transform (in the documentation this is indicated by "If norm='ortho', x[0] and x[N-1] are multiplied by a scaling factor of sqrt(2)") to make sure that Parseval's theorem still holds afterwards. However, then the power in the different modes do not correspond any more to the original signal.

Solution

This modification of the original signal is confusing and I propose the following to anybody who also wants Parseval's theorem to hold while still knowing in which modes the original input signal has power:

f_DCT = scipy.fft.dct(f, type = 1, norm = "backward")

# Apply manual normalisation similar to: norm = "ortho"
# See documentation of SciPy DCT (y[k]).
scaling_factors = np.zeros(np.shape(f_DCT))
scaling_factors[1:-1] = 0.5 * (np.sqrt(2) / np.sqrt(len(f_DCT) -1))
scaling_factors[0] = 0.5 * (1 / np.sqrt(len(f_DCT) -1))
scaling_factors[-1] = 0.5 * (1 / np.sqrt(len(f_DCT) -1))

f_DCT = f_DCT * scaling_factors
del scaling_factors
# Now f_DCT is scaled as expected for norm = "ortho"

# To check Parseval's theorem, one must scale the weight of the first 
# and last data point because of the specific type I of the DCT. 
# See documentation of SciPy DCT (x[0], x[N-1]).
scaling_factors = np.ones(np.shape(f))
scaling_factors[0] = 1 / np.sqrt(2)
scaling_factors[-1] = 1 / np.sqrt(2)

# Compute the signal weighted properly for checking Parseval's theorem (PT).
f_PT = f * scaling_factors
del scaling_factors

# Note that there is now only energy in one single mode (at k=1pi):
array([ 2.93737402e-16,  1.87082869e+00, -2.31984713e-17,  1.78912008e-16,
       -1.78912008e-16,  2.31984713e-17,  0.00000000e+00, -1.25887458e-16])

# Also, Parseval's theorem holds:
np.sum(f_PT * f_PT)   # 3.499999999999999
np.sum(f_DCT * f_DCT) # 3.499999999999999
phil
  • 11
  • 2