dft of sampled sine using python

Question

I'm trying to write a python script to perform a 100-point DFT over a finite length sample of a sinewave at 1/8 the sampling frequency. I'm curious why my DFT magnitude plot has two spikes in it when i plot it.

(I'm wondering if k=0 corresponds to $\omega=0$, k=100 corresponds to $\omega=2\pi$, and k=50 corresponds to $\omega=\pi$... and i'm just seeing an alias imagine above k=50.. or if i'm doing something wrong with DFT calculation..)

# Python 3.5 script
import sys
import numpy as np
import matplotlib.pyplot as plt

pi = np.pi

def stem_xy(x, y, title, filename):
    print("=> STEM_XY")
    fig = plt.figure()
    plt.title(title)
    plt.stem(x, y)
    print("creating file: ", filename)
    fig.savefig(filename)
    plt.close()

def plot_xy(x, y, title, filename):
    print("=> PLOT_XY")
    fig = plt.figure()
    plt.title(title)
    plt.plot(x,y)
    print("creating file: ", filename)
    fig.savefig(filename)
    plt.close()

def dsp_dft_magnitude(DIN, N, filename):
    print("=> DSP_DFT_MAGNITUDE")
    #N: Number of bins
    w_0 = 2*np.pi / N #Fundamental Frequency (Radians)  
    X   = DIN
    XK  = np.zeros(N, dtype=np.complex) 
    XM  = np.zeros(N)

    print("N  = ", N)
    print("W0 = ", w_0)

    for k in range(0,N):
        XK[k]=0
        for n in range(0,len(X)):
            XK[k] = XK[k] + X[n]*np.exp(-1j*k*w_0*n)
        XK[k] = (1/N)*XK[k]
        XM[k] = np.sqrt((XK[k].real)**2 + (XK[k].imag)**2)

    k = range(0,len(XM))
    plot_xy(k, XM, "DFT", filename)

    return XM

def run():
    print("=> RUN")

    sampleF   = 1E3
    sampleN   = 100
    Fsin      = sampleF/8

    sampleT   = 1 / sampleF
    t         = np.linspace(0.0, sampleN*sampleT, sampleN)
    y         = np.sin(2*pi*Fsin*t)

    stem_xy(t, y, "samples", "plot_test0.png")
    dsp_dft_magnitude(y, 100, "plot_test1.png")
    sys.exit()

run()

Answer 1

Found the answer in numpy documents for fft:

# python to perform dft
# from import numpy.fft import *

A = fft(a, n)

A[0] contains the zero-frequency term (the sum of the signal), which is always purely real for real inputs.

A[1:n/2] contains the positive-frequency terms

A[n/2+1:] contains the negative-frequency terms, in order of decreasingly negative frequency.

For an even number of input points, A[n/2] represents both positive and negative Nyquist frequency, and is also purely real for real input.

For an odd number of input points, A[(n-1)/2] contains the largest positive frequency, while A[(n+1)/2] contains the largest negative frequency.

---

There are many ways to define the DFT, varying in the sign of the exponent, normalization, etc. In this implementation, the DFT is defined as

$$ A_k = \sum_{m=0}^{n-1} a_m exp \left\{ -2 \pi i \frac{mk}{n} \right\},\ \ \ \ k=0,...,n-1 $$

The DFT is in general defined for complex inputs and outputs, and a single-frequency component at linear frequency is represented by a complex exponential , where is the sampling interval.

The values in the result follow so-called “standard” order:

Correction to code above:

# python DFT Sin Test
# Bill Moore
# 4/11/2019

import sys
import numpy as np
import matplotlib.pyplot as plt

def dsp_mag_and_phase(DIN, N, filename):
    """
    Cleanup and plots DFT and FFT results

    N        : Number of frequency bins
    filename : png image file to save dft magnitute plot
    """

    print("=> DSP_FREQPLOT")

    XK = DIN
    K  = range(0, N)

    # Swap Positive and Negative Frequencies of DFT/FFT
    #  (to get standard spectrum of -pi to pi)

    # Assign a positive or negative frequency to each point
    KK = []
    for k in K:
        if k>=0 and k <= int(N/2):
            KK.append(k)
        else:
            KK.append(k - int(N))

    #print("KK:", KK)
    #print("len(KK)=", len(KK))
    #sys.exit()

    AM = np.zeros(N)
    AP = np.zeros(N)
    AM = np.zeros(N)
    AP = np.zeros(N)

    #Positive Frequences
    AX_pos = XK[0:int(N/2)+1]
    AK_pos = KK[0:int(N/2)+1]

    #Negative Frequences
    AK_neg = KK[int(N/2)+1:]
    AX_neg = XK[int(N/2)+1:]

    # Concatenate negative then positive 
    AX = np.concatenate((AX_neg, AX_pos))
    AK = np.concatenate((AK_neg, AK_pos))

    # Convert Complex Frequency into Magnitude and Phase
    for k in range(0, len(AK)):
        AM[k] = np.sqrt((AX[k].real)**2 + (AX[k].imag)**2)
        AP[k] = np.arctan2(AX[k].imag, AX[k].real)

    return (AK, AM, AP)


def dsp_dft(DIN, N, filename):
    """
    plots N-point DFT

    N        : Number of frequency bins
    filename : png image file to save dft magnitute plot
    """

    print("=> DSP_DFT_MAGNITUDE")

    DFT = np.zeros(N, dtype=np.complex) 
    K   = range(0, N)
    w0  = 2 * np.pi / N #Fundamental Frequency (Radians)

    # DFT
    for k in K:
        DFT[k]=0
        for n in range(0,len(DIN)):
            DFT[k] = DFT[k] + DIN[n]*np.exp(-1j*k*w0*n)
        DFT[k] = (1/N)*DFT[k]

    # Covert Complex to Magnitude and Phase
    (AK, AM, AP) = dsp_mag_and_phase(DFT, N, filename)

    plot_mag_and_phase(AK, AM, AP, "DFT", filename + ".png")

    return (AK, AM, AP)


def dsp_fft(DIN, N, filename):
    """
    plots N-point FFT 
    (FFT is just a computationally faster DFT)

    N        : Number of frequency bins
    filename : png image file to save dft magnitute plot
    """

    print("=> DSP_FFT_MAGNITUDE")

    # FFT Points equal to length of DIN (default for fft function)
    if N == 0:
        N = len(DIN)

    FFT = np.fft.fft(DIN, N)    

    (AK, AM, AP) = dsp_mag_and_phase(FFT, N, filename)

    plot_mag_and_phase(AK, AM, AP, "FFT", filename + ".png")

    return (AK, AM, AP)

def plot_stem(x, y, title, filename):
    print("=> PLOT_STEM")
    fig = plt.figure()
    plt.title(title)
    plt.stem(x, y)
    print("creating file: ", filename)
    fig.savefig(filename)
    plt.close()

def plot_mag_and_phase(k, mag, phase, title, filename):
    print("=> plot_mag_and_phase")

    fig = plt.figure()
    plt.clf()

    plt.subplot("211")
    plt.title(title)
    plt.ylabel("Magnitude")
    plt.plot(k,mag)

    plt.subplot("212")
    plt.ylabel("Phase")
    plt.plot(k,phase)

    print("creating file: ", filename)
    fig.savefig(filename)
    plt.close()


def dsp_range(sampleN, sampleF):
    """
    return  : time points for sampling at a given frequency

    sampleN : Number of Samples
    sampleF : Sample Frequency in Hz
    """
    sampleT       = 1 / sampleF
    time_range    = np.linspace(0.0, sampleN*sampleT, sampleN)
    return time_range;

def run():
    print("=> RUN")

    sampleF   = 1E3
    sampleN   = 100
    Fsin      = sampleF/4

    # sample function 
    tvec       = dsp_range(sampleN, sampleF)
    y1         = np.sin(2*np.pi*Fsin*tvec)
    plot_stem(tvec, y1, "samples", "plot_test0.png")

    dsp_dft(y1, 100, "plot_test1")
    dsp_fft(y1, 0,   "plot_test2")

    sys.exit()

run()

Answer 2

What you see is negative and positive frequencies. $k=0$ is in the middle, the spikes are at -Fsin and +Fsin.

Answer 3

Its possible to calculate the DFT frequency bins in the correct order from negative frequencies to positive frequencies without needing to rearranging the DFT output at the end.

(In this case the output range corresponding to CFT frequency range $-f_s/2$ to $f_s/2$, or DTFT frequency range $\omega=-\pi$ to $\omega=\pi$.)

N-point DFT:

$$X[\![k]\!] = \sum_{n=-(N/2)+1}^{N/2} x[\![n]\!]\ e^{j2\pi n k / N}$$

When n is less than zero: $$x[\![n]\!] = x[n + N]$$

When n is greater-than or equal to zero: $$x[\![n]\!] = x[n]$$