8
$\begingroup$

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. enter image description here

"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?

$\endgroup$
  • $\begingroup$ To get a deeper understanding of the meaning of a DFT, I recommend you read my first two blog articles: "The Exponential Nature of the Complex Unit Circle" (dsprelated.com/showarticle/754.php) and "DFT Graphical Interpretation: Centroids of Weighted Roots of Unity" (dsprelated.com/showarticle/768.php). $\endgroup$ – Cedron Dawg Apr 27 '18 at 13:57
  • $\begingroup$ Thanks I'll take a look. I'm honestly very surprised at the attention this got when it's all due to a very silly bug in my code. $\endgroup$ – Nimitz14 Apr 27 '18 at 18:20
  • $\begingroup$ I'm surprised too. The continuous vs discrete thing is a big deal though. My blog is all about the discrete case without any reference to the continuous case which is different than teaching the discrete case as a sampled version of the continuous case. $\endgroup$ – Cedron Dawg Apr 27 '18 at 18:44
16
$\begingroup$

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()

enter image description here

$\endgroup$
  • $\begingroup$ btw you were probably having problems because my code assumed python3 ;) $\endgroup$ – Nimitz14 Apr 27 '18 at 18:33
  • 1
    $\begingroup$ @Nimitz14, Not a big deal. I added the "float()" and a bunch of ".0"s on the numbers. Your code ran fine, the only thing I had to remove was the "plt.style.use('ggplot')" statement. $\endgroup$ – Cedron Dawg Apr 27 '18 at 18:46
  • 1
    $\begingroup$ @Nimitz14, I forgot to mention, I added a ft5 routine to the code that pre-calculates the roots of unity values and really shows how the DFT is calculated using the same roots for each bin. $\endgroup$ – Cedron Dawg Apr 27 '18 at 21:19
4
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ If we used i instead of j, we couldn’t say ELI the ICE man. ELJ the JCE man doesn’t have the same ring. Civilization would be threatened $\endgroup$ – Stanley Pawlukiewicz Apr 26 '18 at 3:30
  • 1
    $\begingroup$ elijah the juice man? $\endgroup$ – robert bristow-johnson Apr 26 '18 at 9:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.