I'm coming from an understanding of the continuous-time Fourier Transform, and the effects of doing a DFT and the inverse DFT are mysterious to me.
I have created a noiseless signal as:
import numpy as np
def f(x):
return x*(x-0.8)*(x+1)
X = np.linspace(-1,1,50)
y = f(X)
Now, if I were to perform a continuous Fourier transform on the function $f$ given above, restricted to $[-1,1]$, I would expect the sum of the first few Fourier basis components to give a reasonable approximation to the function $f$ (this is an observation specific to our $f$, since it is approximately sine-wavey over $[-1,1]$). The discrete Fourier transform is an approximation to the continuous one, so assuming that my points y
are sampled noiselessly from $f$ (which they are by design), then the DFT coefficients should approximate the CFT coefficients (I think). So, I obtain a DFT like so (formulae employed):
def DFT(y):
# the various frequencies
terms = np.tile(np.arange(y.shape[0]), (y.shape[0],1))
# the various frequencies cross the equi-spaced "X" values
terms = np.einsum('i,ij->ij',np.arange(y.shape[0]),terms)
# the "inside" of the sum in the DFT formula
terms = y * np.exp(-1j*2*np.pi*terms/y.shape[0])
# sum up over all points in y
return np.sum(terms, axis=1)
def iDFT_componentwise(fy, X):
# this function returns the various basis function components of y, sampled at X
# so the result is a len(X) x len(fy) matrix with each:
# row corresponding to a point in X and each
# column corresponding to a particular frequency.
terms = np.tile(np.arange(len(fy)), (X.shape[0],1))
terms = fy * np.exp(1j*2*np.pi*np.einsum('i,ij->ij',np.arange(X.shape[0])*fy.shape[0]/X.shape[0],terms)/fy.shape[0])
return terms/fy.shape[0]
def iDFT(fy,X):
# summing the Fourier components over all frequencies gives back the original function
return np.sum(iDFT_componentwise(fy,X), axis=1)
I am interested in inspecting the various basis functions that comprise my signal, so I oversample the domain to get a better-resolved picture:
oversampled_X = np.linspace(-1,1,100)
and proceed to check out my components:
fy = DFT(y)
y_f_components = iDFT_componentwise(fy, oversampled_X)
The positive-frequency components look as expected.
import matplotlib.pyplot as plt
plt.plot(oversampled_X, y_f_components[:,1],c='r')
plt.plot(X,y)
plt.show()
However, the negative frequency components look all weird:
plt.plot(oversampled_X, y_f_components[:,49],c='r')
plt.plot(X,y)
plt.show()
This last image looks like it has problems with aliasing. This, in turn, causes problems when I try to reconstitute the function from the Fourier components (see image below)
plt.plot(oversampled_X, iDFT(fy,oversampled_X),c='r')
plt.plot(X,y)
plt.show()
This problem does not occur when I truncate the continuous time Fourier transform of the function to include the same number of terms (see image below):
import sympy
from sympy import fourier_series
from sympy.abc import x
from sympy.utilities.lambdify import lambdify
f = x*(x-0.8)*(x+1)
fourier_f = fourier_series(f, (x, -1, 1))
lambda_fourier_f = lambdify(x,fourier_f.truncate(25),'numpy')
reconstructed_y = lambda_fourier_f(oversampled_X)
plt.plot(oversampled_X,reconstructed_y,c='r')
plt.plot(X,y)
tl;dr
My oversampled inverse Discrete Fourier Transform has a terrible aliasing problem as illustrated here:
The oversampled inverse Discrete Transform:
As opposed to the oversampled inverse Continuous Transform (trucated to the number of terms in the discrete version).
What is the intrinsic property of the DFT that causes this? If the DFT coefficients approximate the CFT coefficients, then why doesn't the CFT have this problem?
Update: The spectrum
As requested, here is the spectrum of $f$. Note that since $f$ is real, the discrete spectrum (excepting the constant term) is symmetric about n/2. I have not attempted to fix the units.
Update2: Extending the function
Per @robertbristow-johnsons suggestion, I decided to check out a slightly different function: $x(x-1)(x+1)$ on $[-1,1]$ (so that the "ends" agree) and I have "repeated" the data a number of times end-to-end. The thought was that this would alleviate some of the weird effects. However, the exact same features appear. (one may wish to open this figure by itself in a new window to enable zooming)
plt.plot(np.abs(fy))
you can see the DFT spectrum. $\endgroup$