SeismologicalDataAnalysis2025/assignment_03/gibbs_incremental_reconstruction.ipynb

In [ ]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from setupFigure import SetupFigure

In [ ]:

def dft_coeff(x):
    """
    Evaluate c_n = 1/N*sum_{k=0}^{N-1}x_k exp(-i*2*pi*nk/N)
    :param x: samples of function with length N
    """
    ns = np.size(x)
    c = np.zeros(ns, dtype=np.complex64)  # zero coefficients before summing
    for n in range(ns):
        for k in range(ns):
            arg = -2*np.pi*1j*n*k/ns
            c[n] = c[n]+x[k]*np.exp(arg)
        c[n] = c[n]/ns
    return c

In the next function, we use the property of the coefficients c_{-n+N}=c_{-n}=c_{n}^*. Thus, we have c_0 and the conjugated pairs c_n, c_{N-n}.

If N is even, then n can run from 1 to N/2-1. The last conjugated pair is c_{N/2-1}, c_{N/2+1}. Finally, we add c_{N/2} which must be real. Thus, we find \begin{align} x_k &= c_0+c_{N/2}\exp(2\pi ik\frac{N/2}{N})+\sum_{n=1}^{N/2-1}[c_n \exp(2\pi i k\frac{n}{N}) + c_{N-n}\exp(2\pi i k\frac{N-n}{N})] \ &= c_0+c_{N/2}\exp(\pi ik)+\sum_{n=1}^{N/2-1}[c_n \exp(2\pi i k\frac{n}{N}) + c_{n}^*\exp(-2\pi i k\frac{n}{N})] \ &= c_0+c_{N/2}\exp(\pi ik)+2\sum_{n=1}^{N/2-1}\mathcal{Re}\left[c_n \exp(2\pi i k\frac{n}{N})\right] ,. \end{align}

For odd N, we again have c_0 and the conjugated pairs c_n, c_{N-n}. Now, n can run from 1 to (N-1)/2. The last conjugated pair is c_{(N-1)/2}, c_{N-(N-1)/2}=c_{(N+1)/2}. There is no coefficient for n=N/2 which should be added. Thus we get: \begin{align} x_k = c_0+2\sum_{n=1}^{(N-1)/2}\mathcal{Re}\left[c_n \exp(2\pi i k\frac{n}{N})\right] ,. \end{align}

In [ ]:

def dft_synthesis(c, nc=0):
    """
    Evaluate the Fourier series as described above
    :param c: complex DFT coefficients
    :param nc: number of conjugated pairs of coefficients to be used 
               (nc <= N/2-1 for even N or nc <= (N-1)/2 for odd N) 
               (default: 0 indicating all coefficients)
    """
    ns = np.size(c)
    
    # initialize x_k with c_0
    x = np.ones(ns)*np.real(c[0])     

    # distinguish even and odd N, we later sum up to nmax
    if ns%2 == 0:                       
        nmax = ns//2-1
    else:
        nmax = (ns-1)//2

    # if N is even and nc == 0, add coefficient c[N/2]*exp(i k pi)
    # the exponential is 1 for even k and -1 for odd k
    if ns%2 == 0 and nc == 0:
        x = x+np.real(c[ns//2])*np.where(np.arange(0, ns)%2 == 0, 1, -1)
    
    # check input nc, reset to nmax if greater
    if nc == 0: nc = nmax
    if nc > nmax: nc = nmax

    # do synthesis
    for n in range(1, nc+1):
        for k in range(ns):
            arg = +2*np.pi*1j*n*k/ns
            x[k] = x[k]+2*np.real(c[n]*np.exp(arg))
    return x

The boxcar function

We now start from a boxcar function which we want to recover from its Fourier coefficients. We do this for an increasing number of coefficients to watch how the series converges.

In [ ]:

tmax = 1000.0
dt = 1.0
ns = int(tmax/dt)
period = ns*dt
t = dt*np.arange(0, ns)
boxcar = np.where(t < 0.3*tmax, 0, 1)*np.where(t > 0.5*tmax, 0, 1)

In [ ]:

fig1, ax1 = SetupFigure(10, 5, "Time [s]", "boxcar", "Boxcar, DT=0.1, TMAX=1000", 14)
ax1.plot(t, boxcar, color='blue', ls='-')

The Fourier coefficients of the boxcar function

We first compute the Fourier coefficients of the boxcar to later carry out an incremental reconstruction.

In [ ]:

c = dft_coeff(boxcar)
f = np.linspace(0, (ns-1)/period, ns)

In [ ]:

nmax = 50
fig2, ax2 = SetupFigure(10, 5, "Frequency [Hz]", "c_n", "Fourier coefficients of Boxcar, DT=1.0, TMAX=1000, Nmax={:d}".format(nmax), 14)
#ax2.plot(f[0:nmax], np.absolute(c[0:nmax]), color='black', ls='-', label='abs')
ax2.plot(f[0:nmax], np.real(c[0:nmax]), color='red', ls='-', label='real')
ax2.plot(f[0:nmax], np.imag(c[0:nmax]), color='blue', ls='-', label='imag')
ax2.legend()

Incremental reconstruction of the boxcar function from the Fourier coefficients

Now we do the Fourier synthesis with an increasing number of coefficients. Watch how the reconstruction changes. What happens at the edges for high numbers of coefficients?

In [ ]:

nc = 20
xr = dft_synthesis(c, nc)
fig3, ax3 = SetupFigure(10, 5, "Time [s]", "xr(t)", "Recovered Boxcar, DT=1, TMAX=1000, Nc={:d}".format(nc), 14)
ax3.plot(t, xr , color='blue', ls='-')
ax3.plot(t, boxcar, color='red', ls=':')

In [ ]: