The following example builds up a square pulse and it applies the functions contained in the module fourier.py to compute the DFT of the function just defined.
Let us import the necessary modules:
import numpy as np
import matplotlib.pyplot as plt
import fourier as f
Now we build the square pulse function, defined as:
\[f(t) = \begin{cases}
0 & \textrm{for}\ 0 \leq t < \frac{\pi}{2} \\
1 & \textrm{for}\ \frac{\pi}{2} \leq t \leq \frac{3\pi}{2} \\
0 & \textrm{for}\ \frac{3\pi}{2} < t. \\
\end{cases}\] Notice that, to avoid the sharp discontinuity at \(t=\pi/2\) and \(t=3\pi/2\) we set the value of the function in these points as the average of the two values \(0\) and \(1\), namely: \(0.5\)!
N=128
dt=2.0*np.pi/N
t=np.arange(0.0, 2.0*np.pi, dt)
y=np.zeros(N)
# Defines the rectangular pulse function in the vector y
for j in range(N):
y[j]=0.0 # Zero everywhere...
if (( t[j] > 0.5*np.pi ) & ( t[j] < 1.5*np.pi )):
y[j]=1.0 # except between pi/2 and 3pi/2...
elif ( t[j] == 0.5*np.pi ):
y[j] = 0.5 # whilst we set 0.5 in pi/2...
elif ( t[j] == 1.5*np.pi ):
y[j] = 0.5 # ... and 3pi/2!
We now plot the vector to be sure the definition was correct:
plt.plot(t,y)
plt.title("A square pulse")
plt.xlabel(r"$t$")
plt.ylabel(r"$y(t)$")
plt.xlim(0.0, 2.0*np.pi)
plt.ylim(-0.2, +1.2)
plt.show()
We now compute and print the Fourier coefficients of the function by using the dft function in the fourier.py module (already included before as "f"!)
cf=f.dft(y,dx)
print(cf)
\[c_n = \begin{cases} 1/2 & \textrm{for}\ n=0 \\ 0 & \textrm{for even}\ n \\ \frac{1}{n\pi} (-1)^{\frac{n-1}{2}} & \textrm{for odd}\ n \\ \end{cases} \]
Now, we compute the spectrum of the function to compare the numerical and theoretical power-law spectrum, by using the function spectrum contained in the module fourier.py:
# Define a vector for the harmonics
nn=np.arange(0,N/2+1)
# Compute the numerical spectrum
sp=f.spectrum(cf,dt)
# Compute the theoretical spectrum
# Define a vector of doubles
thsp = np.ndarray(N/2+1, dtype=np.float64)
for n in nn:
if ( n == 0 ):
thsp[n] = np.pi / 2.0
elif ( np.mod(n,2) == 0 ):
thsp[n] = 0.0
else:
thsp[n] = (4.0/np.pi)*(n**-2.0)
plt.xscale("log")
plt.yscale("log")
plt.xlim(1.0,N/2+1.0)
plt.ylim(-10,1.0)
plt.plot(nn,sp,'b-',thsp,'ro')
plt.show()
Now, we try to reconstruct the function by usign the idft function for several values of the maximum harmonics Nmax:
Nmax=3
ff=f.idft(cf,dt,Nmax)
plt.title("Original and reconstructed function for Nmax="+str(Nmax))
plt.xlabel(r"$t$")
plt.ylabel(r"$y(t)$")
plt.xlim(0.0,2.0*np.pi)
plt.ylim(-0.2,+1.2)
plt.plot(t,y,'b-',t,ff,'ro')
plt.show()
Nmax=11
ff=f.idft(cf,dt,Nmax)
plt.title("Original and reconstructed function for Nmax="+str(Nmax))
plt.xlabel(r"$t$")
plt.ylabel(r"$y(t)$")
plt.xlim(0.0,2.0*np.pi)
plt.ylim(-0.2,+1.2)
plt.plot(t,y,'b-',t,ff,'ro')
plt.show()
Nmax=41
ff=f.idft(cf,dt,Nmax)
plt.title("Original and reconstructed function for Nmax="+str(Nmax))
plt.xlabel(r"$t$")
plt.ylabel(r"$y(t)$")
plt.xlim(0.0,2.0*np.pi)
plt.ylim(-0.2,+1.2)
plt.plot(t,y,'b-',t,ff,'ro')
plt.show()
Nmax=64
ff=f.idft(cf,dt,Nmax)
plt.title("Original and reconstructed function for Nmax="+str(Nmax))
plt.xlabel(r"$t$")
plt.ylabel(r"$y(t)$")
plt.xlim(0.0,2.0*np.pi)
plt.ylim(-0.2,+1.2)
plt.plot(t,y,'b-',t,ff,'ro')
plt.show()
rfft(A), irfft(A), rfft2(A), irfft2(A), rfftn(A), irfftn(A)which compute the FFT of a real 1, 2, n-dimensional numpy.ndarray-s.
fft(A), ifft(A), fft2(A), ifft2(A), fftn(A), ifftn(A)We are interested, for the moment, to 1-dimensional transforms, therefore to the functions rfft(A) and irfft(A). We notice that:
Example: let us consider the sawtooth function \(s(t)\), defined as: \[ s(t) =
\begin{cases}
t & \textrm{for}\ 0 \leq t < \pi \\
0 & \textrm{for}\ t = \pi \\
t - 2\pi & \textrm{for}\ \pi < t \leq 2\pi \\
\end{cases}
\] and let us apply the numpy.fft functions to compute the Fourier coefficients and to reconstruct the function:
# Define the sawtooth function s(t)
N2 = int( N / 2 )
s=np.zeros(N, dtype=np.float64)
s[0:N2] = t[0:N2]
s[N2] = 0.0
s[N2+1:] = t[N2+1:] - 2.0*np.pi
plt.title("Sawtooth function")
plt.xlabel(r"$t$")
plt.ylabel(r"$s(t)$")
plt.plot(t,s,'-')
plt.show()
# Computes the Fourier coefficients of the function through the python rfft function
cfs = np.fft.rfft(s)
# Normalize the coefficients by dividing for N and print them...
cfs = cfs / N
print(cfs)
We can now use the function spectrum from the fourier.py module, to calculate the spectrum of the function and to compare with a spectrum \(n^{-2}\) again:
# Compute the numerical spectrum
sps=f.spectrum(cfs,dt)
# Plot the spectrum along with n^-2 function
plt.title(r"Spectrum of the sawtooth function $s(t)$")
plt.xlabel(r"$n$")
plt.ylabel(r"$E(n)$")
plt.xscale("log")
plt.yscale("log")
plt.xlim(1.0,N/2+1.0)
plt.ylim(1.e-10,100.0)
# Notice that we multiply n^-2 for a factor 10, to
plt.plot(nn,sps,'b-',10.0*nn**(-2.0),'ro')
plt.show()
Here is the function reconstruct:
def reconstruct( ak, n ):
Nyq = ak.size - 1 # Nyquist's number
Ns = 2 * Nyq # Number of points in the output vector
sr = np.zeros( Ns, dtype=np.float64 ) # Initialize the output vector with zeros
if ( n > Nyq ):
print("Error! The number of harmonics for the reconstruction must be lesser than the Nyquist number!")
return sr
else:
akk = np.zeros(Nyq+1, dtype=np.complex128) # Define a new vector for the Fourier's coefficients
akk[:n+1] = Ns * ak[:n+1] # Copy the Fourier's coefficients up to n, into the new vector
# Note the product for Ns to re-normalize the coefficients
sr = np.fft.irfft(akk) # Inverse Fourier transform of the reconstructed function
return sr
Nmax=5
plt.title("Sawtooth function and reconstructed function for Nmax = " + str(Nmax))
plt.xlabel(r"$t$")
plt.ylabel(r"$s(t), s(t,N_{max})$")
plt.xlim(0.0,2.0*np.pi)
plt.ylim(-np.pi,+np.pi)
srec=reconstruct(cfs,Nmax)
plt.plot(t,s,'b-',t,srec,'ro')
plt.show()
Nmax=11
plt.title("Sawtooth function and reconstructed function for Nmax = " + str(Nmax))
plt.xlabel(r"$t$")
plt.ylabel(r"$s(t), s(t,N_{max})$")
plt.xlim(0.0,2.0*np.pi)
plt.ylim(-np.pi,+np.pi)
srec=reconstruct(cfs,Nmax)
plt.plot(t,s,'b-',t,srec,'ro')
plt.show()
Nmax=32
plt.title("Sawtooth function and reconstructed function for Nmax = " + str(Nmax))
plt.xlabel(r"$t$")
plt.ylabel(r"$s(t), s(t,N_{max})$")
plt.xlim(0.0,2.0*np.pi)
plt.ylim(-np.pi,+np.pi)
srec=reconstruct(cfs,Nmax)
plt.plot(t,s,'b-',t,srec,'ro')
plt.show()
Nmax=64
plt.title("Sawtooth function and reconstructed function for Nmax = " + str(Nmax))
plt.xlabel(r"$t$")
plt.ylabel(r"$s(t), s(t,N_{max})$")
plt.xlim(0.0,2.0*np.pi)
plt.ylim(-np.pi,+np.pi)
srec=reconstruct(cfs,Nmax)
plt.plot(t,s,'b-',t,srec,'ro')
plt.show()