1
$\begingroup$

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()

enter image description here

However, the negative frequency components look all weird:

plt.plot(oversampled_X, y_f_components[:,49],c='r')
plt.plot(X,y)
plt.show()

enter image description here

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:

enter image description here

As opposed to the oversampled inverse Continuous Transform (trucated to the number of terms in the discrete version).

enter image description here

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.

enter image description here

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)

enter image description here

$\endgroup$
25
  • 3
    $\begingroup$ first time i saw someone put "tl;dr" in their post!! :-) $\endgroup$ Aug 26, 2019 at 21:39
  • 3
    $\begingroup$ //"What is the intrinsic property of the DFT that causes this?"// --- what is "this"? the glitches at the two sides? $\endgroup$ Aug 26, 2019 at 21:49
  • 2
    $\begingroup$ Are you trying to interpolate in the time domain (I.e same time duration, higher sample rate) or in the frequency domain (longer time duration, same sample rate)? Also, you are doing your interpolation and inverse DFT in one step. If you could show what your frequency data looks like before and after the interpolation, it would be easier to help. $\endgroup$
    – Dan Szabo
    Aug 26, 2019 at 23:34
  • $\begingroup$ interpolating in the time domain. "Doing interpolation and inverse DFT in one step" yes, because I'm simply mirroring the inverse CFT process. The inverse CFT gives values "near" the curve so long as I am interpolating, no matter where I evaluate it. The question is why doesn't the inverse DFT have this same property? I'll post the frequency-domain pics when I get a chance. In the meantime, if you plt.plot(np.abs(fy)) you can see the DFT spectrum. $\endgroup$
    – Him
    Aug 27, 2019 at 3:05
  • 3
    $\begingroup$ it looks like you're using Python. sorry, but i haven't learned Python, so i cannot examine your code. maybe @Fat32 knows what's going on. my suspicion is that perhaps you are not treating the top half of the DFT output as "negative frequency" components. i dunno. but somehow an image at twice the original sample rate is getting in there and making it look like you upsampled by a factor of two by zero inserting between each original sample. if you're doing that mistake, your result should not be purely real, but should have a significant imaginary component. $\endgroup$ Aug 27, 2019 at 20:04

2 Answers 2

1
$\begingroup$

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)

enter image description here

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.

$\endgroup$
6
  • $\begingroup$ it's a miracle that you could infer Eq.2 from this post. Nowhere it ever comes close to imply that... ;-) rather it loosely talks about densely sampling on some unknown domain grid, resulting in weird positive and negative frequency components and aliasing somewhere in some domains(?) $\endgroup$
    – Fat32
    Aug 28, 2019 at 9:16
  • $\begingroup$ @Fat32: That's basically what the OP's function iDFT_componentwise() does, as far as I can see. $\endgroup$
    – Matt L.
    Aug 28, 2019 at 11:43
  • 1
    $\begingroup$ Equation 2 is exactly what I'm doing. I had interpreted it, however, as plugging $n \in [0,N-1]$ into Equation (1), but at non-integral points (i.e. instead of $n \in Z_{N-1}$ I was plugging in $n$ at arbitrary points $\in [0,N-1]$). Equation (2) shuffles the fraction out of $n$ and into the normalization factor ($/N$) to recast the problem in a different light. "Note that these values are not even real-valued" I had not previously noticed this, but this is true. I will mull all of this over today. $\endgroup$
    – Him
    Aug 28, 2019 at 13:19
  • $\begingroup$ AHHHHHHHH..... Since my original data have finite support, the Fourier transform is actually periodic! If I want, instead, to pretend that my original data were periodic, I need to interpret the computed Fourier coefficients a being a single period of the Fourier transform. The single period that is centered at zero corresponds to exactly this imaginary scenario of my original data having been a single period of some continued function. Cool stuff. $\endgroup$
    – Him
    Aug 28, 2019 at 14:48
  • 1
    $\begingroup$ Ah. It is because truncating the otherwise periodic Fourier transform to that particular period is actually the Fourier transform of some complex function or another that just happens to agree with my $f$ on the original sampled datapoints. Thanks. This answer was invaluable in bringing me to a comprehension of the matter. $\endgroup$
    – Him
    Aug 28, 2019 at 15:32
0
$\begingroup$

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 :

enter image description here

$\endgroup$
12
  • $\begingroup$ "oversampled CFT approximation." It is not an approximation. Since I know the function from which I sampled the data, I can perform a symbolic Fourier transform, and evaluate that at any desired points on the support. $\endgroup$
    – Him
    Aug 27, 2019 at 19:39
  • $\begingroup$ ok whatever. This is my answer to the image you pointed in your comment. There are zero samples for every other and the answer is as above, I don't know what you are doing wrong, furthermore, you also say that you insert those zeros which are actually nonzero DC constants ? that's quite strange and needs further explanation. $\endgroup$
    – Fat32
    Aug 27, 2019 at 19:50
  • $\begingroup$ I am not inserting any zeros. The inverse discrete Fourier transform is a function of time. I am simply plugging an evenly-sampled grid on [-1,1] of 100 points into the inverse discrete Fourier transform. Zero is what's coming out. $\endgroup$
    – Him
    Aug 27, 2019 at 19:52
  • $\begingroup$ (they're actually not zero, they're a constant value, the first DFT coefficient). are they zero or nonzero ? $\endgroup$
    – Fat32
    Aug 27, 2019 at 19:53
  • $\begingroup$ *a constant is what's coming out. $\endgroup$
    – Him
    Aug 27, 2019 at 19:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.