CODEWITHSUNDEEP
pythonfftcomplex-numbers
Decimak
Decimak

Reputation: 35

Discrete fourier transform giving complex conjugate of "right" answer

I'm attempting to implement a simple version of the discrete fourier transform in python. My code is as follows: 

#!/usr/bin/env python
import cmath
def dft_simple(sequence):
# dft of seq defined as
# sigma from n=0 to N-1 of x(n) *exp(-2*pi*j*k*n/N)
  seqLenth = len(sequence)
  complexSequence = []
  for k in range(seqLenth):
    sigma = 0 - 0j
    print("k = {}".format(k))
    for n in range(seqLenth):
      print("n = {}".format(n))
      print("value = {}".format(sequence[n] * cmath.exp(-2*1j * cmath.pi * float(k) \
                                    * float(n)  / float(seqLenth))))
      sigma = sigma + (sequence[n] * cmath.exp(-2*1j * cmath.pi * float(k) \
                                    * float(n)  / float(seqLenth)))
      print("exp = {0}".format(-2*1j * cmath.pi * float(k) \
                                    * float(n)  / float(seqLenth)))
    complexSequence.append(sigma)

    print("sum = {}".format(sigma))
    print("")
  return(complexSequence)
seq4 = [1,1,1,1,0,0,0,0]
print(dft_simple(seq4))

I receive a result of: 

[(4+0j), (1-2.414213562373095j), (-1.8369701987210297e-16-2.220446049250313e-16j), (1-0.4142135623730949j), -2.449293598294706e-16j, (0.9999999999999992+0.4142135623730959j), (3.2904645469127765e-16-3.3306690738754696e-16j), (0.9999999999999997+2.4142135623730954j)]

This differs from the answer I get on wolfram alpha when calculating the DFT of the same sequence here, in two ways. First, wolfram alpha is dividing by sqrt(N), where N is the length of the sequence, which is just a different symmetric definition of the forward and inverse transforms.
Secondly, and more confusingly, my implementation is giving me the complex conjugate of the result wolfram alpha is giving me - the numeric values are otherwise approximately the same. Is this an issue of an implementation problem on my code(e.g. syntax error), or is it simply wolfram using a different definition of the discrete fourier transform?

Upvotes: 3

Views: 1045

Answers (1)

SleuthEye
SleuthEye

Reputation: 14579

In both cases (for the scaling and for the complex conjugated result), the underlying cause is the difference in the definition used for the Discrete Fourier Transform (DFT).

The default definition of the DFT from Wolfram uses the formula:

\frac{1}{n^{1/2}}\sum_{r=1}^n u_r e^{2\pi i (r-1)(s-1)/n}

Or equivalently using zero based indexing, time index n, frequency index k and j=sqrt(-1) to compare against your implementation:

\frac{1}{N^{1/2}}\sum_{n=0}^{N-1} u_n e^{2\pi j k n/N}

Your implementation uses what Wolfram refers to as the "signal processing" convention:

\sum_{r=1}^n u_r e^{-2\pi i (r-1)(s-1)/n}

which again is equivalent to:

\sum_{n=0}^{N-1} u_n e^{-2\pi j k n/N}

For real-valued input sequences, using a negative sign in the complex exponential term would produce a result that is the complex conjugate of the similar expression which uses a positive sign in the complex exponential term (and vice versa):

\begin{align}\sum_{n=0}^{N-1} u_n e^{-2\pi j k n/N}&= \sum_{n=0}^{N-1} u_n \mbox{conjugate}(e^{2\pi j k n/N}) \ &= \sum_{n=0}^{N-1} \mbox{conjugate}(u_n e^{2\pi j k n/N}) &,\mbox{for real $u_n$}\ &= \mbox{conjugate}\left(\sum_{n=0}^{N-1} u_n e^{2\pi j k n/N}\right) &,\mbox{for real $u_n$}\\end{align}

Upvotes: 2

Related Questions