Aliasing when interpolating with DFT?

Question

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)

Answer 1

Let me summarize my understanding of what you're trying to do. You have a real-valued sequence $x[n]$, obtained by sampling a real-valued continuous function, and you computed its DFT $X[k]$. The sequence can be expressed in terms of its DFT coefficients:

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N},\qquad n\in[0,N-1]\tag{1}$$

where $N$ is the length of the sequence.

Now you want to interpolate that sequence, and I believe you're trying to do this in the following way:

$$\tilde{x}[m]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi mk/M},\qquad m\in[0,M-1],\quad M>N\tag{2}$$

This, however, doesn't work. If $M$ happens to be an integer multiple of $N$, then $\tilde{x}[nM/N]=x[n]$ is satisfied, but the other values of $\tilde{x}[m]$ are by no means interpolated values of $x[n]$. Note that these values are not even real-valued.

What you can do is approximately compute the Fourier coefficients of the original continuous function using the (length $N$) DFT of the sampled function, and then approximately reconstruct samples of the function on a dense grid (of length $M>N$):

$$\tilde{x}[m]=\frac{1}{N}\sum_{k=-K}^KX[k]e^{j2\pi mk/M},\qquad m\in[0,M-1]\tag{3}$$

Note that in $(3)$ the summation indices are symmetric, and the number $K$ cannot exceed $N/2$ because that's the number of independent DFT coefficients you have due to conjugate symmetry of $X[k]$ (because $x[n]$ is assumed to be real-valued).

Eq. $(3)$ is just equivalent to zero-padding in the frequency domain, which corresponds to interpolation in the time domain. Note, however, that the zero padding is done in such a way that conjugate symmetry is retained, i.e., the zeros are inserted around the Nyquist frequency, and not simply appended to the DFT coefficients.

With $X[-k]=X[N-k]$ and $X[k]=X^*[N-k]$, Eq. $(3)$ can be rewritten as

$$\begin{align}\tilde{x}[m]&=\frac{1}{N}X[0]+\frac{1}{N}\sum_{k=1}^K\left(X[k]e^{j2\pi mk/M}+X[-k]e^{-j2\pi mk/M}\right)\\&=\frac{1}{N}X[0]+\frac{1}{N}\sum_{k=1}^K\left(X[k]e^{j2\pi mk/M}+X^*[k]e^{-j2\pi mk/M}\right)\\&=\frac{1}{N}X[0]+\frac{2}{N}\textrm{Re}\left\{\sum_{k=1}^KX[k]e^{j2\pi mk/M}\right\},\qquad m\in[0,M-1]\end{align}\tag{4}$$

The following Matlab/Octave code illustrates the above:

N = 100;
t = linspace (-1,1,N);

M = 200;
ti = linspace (-1,1,M);

x = t .* (t - 0.8) .* (t + 1);
x = x(:);

X = fft(x);
X = X(:);

Nc = 20;    % # Fourier coefficients (must not exceed N/2)

x2 = X(1) + 2*real( exp( 1i * 2*pi/M * (0:M-1)' * (1:Nc-1) ) * X(2:Nc) );
x2 = x2 / N;

plot(t,x,ti,x2)

Note that the approximation of the blue curve by the green curve in the above figure is two-fold: first, there's only a finite number of Fourier coefficients, and second, the Fourier coefficients are only approximately computed from samples of the original function.

Answer 2

I have looked at your image (and didn't read the post) and its explanation is as follows.

Let the data in the blue curve be $x[n]$, and the data in the red curve be $y[n]$; then it can be seen and shown that:

$$y[n] = \tfrac12 ( x[n] + (-1)^n x[n] ) $$

In the DTFT domain this relationship becomes :

$$ Y(e^{j\omega}) = \tfrac12 \big( X(e^{j\omega}) + X(e^{j(\omega-\pi)}) \big) $$

In the DFT domain, for even $N$, this becomes: $$ Y[k] = \tfrac12 \big( X[k] + X[k - \tfrac{N}2] \big) $$

Now, in the frequency domain we do not have an interpolated spectrum but an aliased one... Who does or what causes this aliasing? I don't know as I didn't read your post in such detail, but it's neverthess there...

Hence when you convert $Y[k]$ back into time-domain, you will obtain $y[n]$ sequence as an inversion of an aliased spectrum of $X[k]$ ; There's no interpolation in any domains...

The following MATLAB / OCTAVE code demonstrates what you try to achieve by your oversampled CFT approximation(?) of DFT.

L = 32;
t = linspace(-1,1,L);
x = t.*(t-0.8).*(t+1);          % Polynomial signal of length L

N = L ;                         % DFT length... (make it even) 
X = fft(x,N);

Y = 0.5*(X + fftshift(X));      % is this an "INTERPOLATION " of
                                % X[k]? No. It's: X[k] + X[k-N/2]

y = real(ifft(Y,N));            % inverse-DFT to reconstruct y[n]
y = y(1:L);                     % Trim the trailing zeros, if any, on the
                                % reconstructed y[n] of length N

figure,stem(t,x)
hold on
stem(t,y,'+r');title('x[n] vs y[n]');   % 
legend('x[n]','','y[n]');

The result is :