DFT and perfect reconstruction of a square wave on a digital computer

Question

I know that in theory, when reconstructing a square wave from its Fourier coefficients, unless we have an infinite amount of them, the resulting reconstruction will have Gibbs ringing artifacts due to lack of enough harmonics.

On a computer, we can take the Fourier transform X = fft(x) of a square wave x, and reconstruct it without artifact with x_rec = ifft(X), maybe with some rounding error of the order of 1e-17 or something but no visible ringing.

I don't have a satisfying answer for that? I guess there has to be something to do with the fact the "the square wave" x is a digitized version of a continuous wave, and my Fourier basis vectors (complex exponentials are also discretized of course since we are in a computer...) but still... how would you justify the absence of Gibbs ringing artefacts from Fourier reconstruction of the Fourier transform of a digital square wave ?

%%%%%%%%%%%%%%%%%%%% Tought experiment proposed by  Dan Szabo
fs=10;%sampling frequency
t=0:(1/fs):1-(1/fs);


s = [1 1 1 1 1 0 0 0 0 0];
sTr = imtranslate(s,[0.5 0])

sTr =

    0.5000    1.0000    1.0000    1.0000    1.0000    0.5000         0         0         0         0

Answer 1

First of all, it's a misunderstanding that the Gibbs phenomenon disappears if you use infinitely many Fourier series coefficients to reconstruct a discontinuous periodic function, such as a square wave. It doesn't. The reason is that generally the Fourier series doesn't converge point-wise, but it converges in the mean, i.e.,

$$\lim_{N\to\infty}\int_{0}^{T}\left|x(t)-\sum_{k=-N}^{N}c_ke^{j2\pi kt/T}\right|^2= 0\tag{1}$$

if $x(t)$ is a $T$-periodic function, and $c_k$ are its Fourier coefficients.

Taking the discrete Fourier transform (DFT) of a finite sequence of numbers just corresponds to a matrix multiplication:

$$\mathbf{y}=\mathbf{Ax}\tag{2}$$

and as long as the matrix $\mathbf{A}$ is invertible you can compute $\mathbf{x}$ from $\mathbf{y}$:

$$\mathbf{x}=\mathbf{A}^{-1}\mathbf{y}\tag{3}$$

This has nothing to do with the Fourier series of a continuous periodic function, and it has nothing to do with the Gibbs phenomenon.

Answer 2

You can certainly create a discrete square wave and take it's FFT. However, the result will NOT be the spectrum of a square wave. Since the square wave is not band limited, you will get substantial amount of aliasing. Keep in mind that everything you represent numerically in a computer is discrete and periodic in both time and frequency.

In order to match the spectrum of the square wave, you would have to low-pass filter your time domain signal first, and that would certainly generate a lot of ringing.

Answer 3

This is easier understood by the principles of Linear Algebra. The DFT is a lossless linear transform. The span of basis vectors is the full vector space, so it doesn't matter what your signal is, you can always reconstruct it losslessly.

---

Response to the last comment.

Leaving the ranges off, each is a dot product corresponding to calculating a single element of your output vector.

$$ X[k] = \sum x[n] e^{-i 2\pi \frac{n}{N} k } $$

$$ F(\Omega) = \int f(t) e^{-i 2\pi t \Omega } $$

In the discrete case, $k$ is the element index corresponding to frequency. In the continuous case $\Omega$ is in the frequency domain. They are both just dot products, aka inner products.

Sweep across the $k$s, you have a matrix multiplication. Sweep across the $\Omega$s, the same, only continuous.

Think of the basis vector as this:

$$ g_{\Omega}(t) = e^{-i 2\pi t \Omega } $$

Then:

$$ F(\Omega) = f \cdot g_{\Omega} $$

I can't make it clearer than that. Nor am I going to try.