Implementation of Wikipedia Equation for the DFT

Question

I was writing a simple fourier transform implementation and looked at the DFT equation on wikipedia for reference, when I noticed that I was doing something differently, and after thinking about it felt that the wikipedia version must be wrong because it's very simple to think of a signal that when fourier transformed (with that equation) will return an incorrect spectrum: Because the equation wraps the signal around the complex plane only once (due to the $n/N$ with $0<n<N-1$), any signal that is periodic an even number of times (while wrapping the complex plane) will have no spectrum as the usual peaks (while going around the unit circle) that would appear during a DFT will cancel each other out (when an even number of them appear).

To check this I wrote some code which produced the following image, which seems to confirm what my thoughts.

"Time using equation" uses the equation $$ X_f = \sum_{n=0}^{N-1} x_n (\cos(2\pi f t_n) - i \sin(2\pi f t_n)) $$ with $t$ a vector of time (so the time $t_n$ at which $x_n$ was sampled for example). It can be found in the function ft below.

The wikipedia equation, linked above, is copied here for reference: $$ X_f = \sum_{n=0}^{N-1} x_n \left(\cos\left(2\pi f \frac{n}{N}\right) - i\sin\left(2\pi f \frac{n}{N}\right)\right) $$ It can be found in the function ft2.

import numpy as np
import matplotlib.pyplot as plt 
plt.style.use('ggplot')

def ft(t, s, fs):
    freq_step = fs / len(s)
    freqs = np.arange(0, fs/2 + freq_step, freq_step)
    S = []
    for freq in freqs:
        real = np.sum(s * np.cos(2*np.pi*freq * t)) 
        compl = np.sum(- s * np.sin(2*np.pi*freq * t)) 
        tmpsum = (real**2 + compl**2) ** 0.5 
        S.append(tmpsum)
    return S, freqs

def ft2(s, fs):  # Using wikipedia equation
    nump=len(s)
    freq_step = fs / nump
    freqs = np.arange(0, fs/2 + freq_step, freq_step)
    S = []
    for i, freq in enumerate(freqs):
        real = np.sum(s * np.cos(2*np.pi*freq * i/nump))
        compl = np.sum(- s * np.sin(2*np.pi*freq * i/nump))
        tmpsum = (real**2 + compl**2) ** 0.5 
        S.append(tmpsum)
    return S, freqs


def main():
    f = 5 
    fs = 100 
    t = np.linspace(0, 2, 200)
    y = np.sin(2*np.pi*f*t) + np.cos(2*np.pi*f*2*t)

    fig = plt.figure()
    ax = fig.add_subplot(311)
    ax.set_title('Signal in time domain')
    ax.set_xlabel('t')
    ax.plot(t, y)

    S, freqs = ft(t, y, fs) 

    ax = fig.add_subplot(312)
    ax.set_xticks(np.arange(0, freqs[-1], 2)) 
    ax.set_title('Time using equation')
    ax.set_xlabel('frequency')
    ax.plot(freqs, S)

    S, freqs = ft2(y, fs) 
    ax = fig.add_subplot(313)
    ax.set_title('Using Wiki equation')
    ax.set_xlabel('frequency')
    ax.set_xticks(np.arange(0, freqs[-1], 2)) 
    ax.plot(freqs, S)

    plt.tight_layout()
    plt.show()

main()

Obviously it seems rather unlikely that I would have randomly found an error on such a high profile wiki page. But I can't see a mistake in what I've done?

Answer 1

You have a bug in ft2. You are incrementing i, and freq together. That's not how you want your summation to work. I messed around with fixing it, but it got messy. I decided to rewrite it from a discrete perspective instead of trying to use the continuous terminology. In the DFT, the sampling rate is irrelevant. What matters is how many samples are used (N). The bin numbers (k) then correspond to frequency in units of cycles per frame. I tried to keep you code as intact as possible so it would remain easily comprehensible to you. I also unfurled the DFT calculation loops to hopefully reveal their nature a little bit better.

Hope this helps.

Ced

import numpy as np
import matplotlib.pyplot as plt 

def ft(t, s, fs):
    freq_step = fs / len(s)
    freqs = np.arange(0, fs/2, freq_step)
    S = []
    for freq in freqs:
        real = np.sum(s * np.cos(2*np.pi*freq * t)) 
        compl = np.sum(- s * np.sin(2*np.pi*freq * t)) 
        tmpsum = (real**2 + compl**2) ** 0.5 
        S.append(tmpsum)
    return S, freqs

def ft3(s, N):  # More efficient form of wikipedia equation

    S = []

    slice  = 0.0
    sliver = 2*np.pi/float(N) 

    for k in range(N/2):

        sum_real = 0.0    
        sum_imag = 0.0
        angle = 0.0
        for n in range(N):
            sum_real +=  s[n] * np.cos(angle)
            sum_imag += -s[n] * np.sin(angle)
            angle += slice

        slice += sliver
        tmpsum = (sum_real**2 + sum_imag**2) ** 0.5 
        S.append(tmpsum)

    return S

def ft4(s, N):  # Using wikipedia equation

    S = []

    for k in range(N/2):

        sum_real = 0.0    
        sum_imag = 0.0
        for n in range(N):
            sum_real +=  s[n] * np.cos(2*np.pi*k*n/float(N))
            sum_imag += -s[n] * np.sin(2*np.pi*k*n/float(N))

        tmpsum = (sum_real**2 + sum_imag**2) ** 0.5 
        S.append(tmpsum)

    return S

def ft5(s, N):  # Roots of Unity Weighted Sum

    sliver = 2 * np.pi / float( N ) 

    root_real = np.zeros( N )
    root_imag = np.zeros( N )

    angle = 0.0
    for r in range(N):
        root_real[r] =  np.cos( angle )
        root_imag[r] = -np.sin( angle )
        angle += sliver

    S = []

    for k in range( N/2 ):

        sum_real = 0.0    
        sum_imag = 0.0
        r = 0

        for n in range( N ):
            sum_real += s[n] * root_real[r]
            sum_imag += s[n] * root_imag[r]
            r += k
            if r >= N : r -= N

        tmpsum = np.sqrt( sum_real*sum_real + sum_imag*sum_imag )
        S.append( tmpsum )

    return S

def main():

    N  = 200
    fs = 100.0

    time_step = 1.0 / fs
    t = np.arange(0, N * time_step, time_step)

    f = 5.0
    y = np.sin(2*np.pi*f*t) + np.cos(2*np.pi*f*2*t)

    fig = plt.figure()
    ax = fig.add_subplot(311)
    ax.set_title('Signal in time domain')
    ax.set_xlabel('t')
    ax.plot(t, y)

    S, freqs = ft(t, y, fs) 

    ax = fig.add_subplot(312)
    ax.set_xticks(np.arange(0, freqs[-1], 2)) 
    ax.set_title('Time using equation')
    ax.set_xlabel('frequency')
    ax.plot(freqs, S)

    S = ft3(y, N) 
    ax = fig.add_subplot(313)
    ax.set_title('Using Wiki equation')
    ax.set_xlabel('frequency')
    ax.set_xticks(np.arange(0, freqs[-1], 2)) 
    print len(S), len(freqs)
    ax.plot(freqs, S)

    plt.tight_layout()
    plt.show()

main()

Answer 2

i am not gonna look through your code. the wikipedia page looks okay, but it is a good example of the "format war" or "notation war" or "style war" between mathematicians and electrical engineers. some of it, i think the math people are right. EEs should have never adopted "$j$" for the imaginary unit. that said, this is a better expression of the DFT and inverse is:

DFT: $$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi nk/N} $$

iDFT: $$ x[n] = \frac{1}{N} \sum\limits_{k=0}^{N-1} X[k] \, e^{j 2 \pi nk/N} $$

because electrical engineers doing DSP like to use $x[n]$ as the a sequence of samples in "time" and $X[k]$ as the sequence of discrete samples in "frequency". mathematicians might like this better:

DFT: $$ X_k = \sum\limits_{n=0}^{N-1} x_n \, e^{-i 2 \pi nk/N} $$

iDFT: $$ x_n = \frac{1}{N} \sum\limits_{k=0}^{N-1} X_k \, e^{i 2 \pi nk/N} $$

and that is the same as the wikipedia page.

you might need to pay more attention to the use of $+$ or $-$ in the exponent and how that translates to $+$ or $-$ against the $\sin(\cdot)$ term.

Answer 3

I came back to this and tried deriving the discrete version which helped make things more sense:

Somehow $f_k t_n = f(n, k, N)$

$ f_k = \frac{f_s}{N}k$ and $ t_n = \frac{T}{N}n $

$f_s = \frac{N}{T}$

So

$f_k t_n = \frac{f_s}{N}k\frac{T}{N}n = \frac{N}{TN}k\frac{T}{N}n = \frac{kn}{N}$

Done!

edit2: I have just realized that the visual reason the result is similar for when k=n and k=N-n is because in one case you're just going around the circle clockwise, and in the other case counterclockwise. So of course the real (horizontal) component will be the same while the imaginary (vertical) component will be inversed.