Implementation of Wikipedia Equation for the DFT 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&lt;n&lt;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?
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()