# dft of sampled sine using python

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()  Found the answer in numpy documents for fft:

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

A = fft(a, n)


A 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() What you see is negative and positive frequencies. $$k=0$$ is in the middle, the spikes are at -Fsin and +Fsin.

• "k=0 to k=n/2" corresponds to "$\omega=0$ to $\omega=\pi$" and "k=n/2+1 to k=n" corresponds to "$\omega=-\pi$ to $\omega=0$"? In this case, I can safely, throw away half of my DFT points without any lost of information... – Bill Moore Apr 10 '19 at 17:17

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]$$