import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
def theOriginalFunction(var,period):
return np.piecewise(var, [((var >= 0)*(var < period/2)), ((var >= period/2)*(var < period))], [1, -1])
def innerSum(j,Yn,N):
inner = np.array([Yn[n]*np.exp(-2*np.pi*1j*j*n/N) for n in range((-N/2)+1,N/2)])
return inner.sum()
#total = 0
#for n in range (1, N):
# total = total + (Yn[n]*np.exp(-2*np.pi*1j*((j*n)/N)))
#return total
# Generate a function with period L
L = 2*np.pi
x = np.linspace(-10, 10, 1000)
y = theOriginalFunction(x,L)
#Generate the corresponding discrete Fourier transform
N = 128
xn = np.linspace(-10, 10, N)
Yn = theOriginalFunction(xn,L)
#myTest = innerSum(j,Yn,N)
#print myTest
#Cj = np.array([innerSum(j,Yn,N).imag for j in range(0, N-1)])
Cj = np.empty(N) + 1j * np.empty(N)
for j in range(0, N-1):
Cj[j] = (L/N)*innerSum(j,Yn,N)
Cj_abs = np.sqrt(Cj.real**2 + Cj.imag**2)
#Generate the plot
plt.plot(x, y, 'k-')
#plt.plot(xn, Cj.real)
#plt.plot(xn, Cj.imag)
plt.plot(xn, Cj_abs)
plt.margins(0.1)
plt.show()
Cj_abs[:10]
Cj_abs[-10:]