$L2$ and BandLimted Function Spaces $V$

$L2$ space is defined as space of square-integrable equivalent classes of functions: \(L2 = \bigg\{f: \|f\|_{L2} = (\int_{-\infty}^{\infty} \|f(t)\|^ dt)^{1/2} dt < \infty \bigg\}\) It is endowed with inner product of the form \(\langle f,g\rangle_{L2} = \int_{-\infty}^{\infty} f(t)\cdot g(t) dt\)

Bandlimited Space $V$ is a closed subspace of $L2$ defined as

\(V = L2 \cap B_{[-\Omega/2, \Omega/2]}\) where functions in $B_{[-\Omega/2, \Omega/2]}$ are characterized by their fourier transform having support only in domain $\omega \in [-\Omega/2, \Omega/2]$. \(B = \bigg\{f: Supp\{\hat{f}(\xi)\} = [-\Omega/2, \Omega/2] \bigg\}\) where \(\hat{f}(\xi) = \int_{\mathbb{R}^{d}} f(t)e^{2\pi i \langle \xi,t\rangle} dt\) is the fourier transform of $f$.

Analysis -> Synthesis as Approximation on Projected Space

Sampling can also be considered as analysis, where measurement (in $\ell2$) of the singal in $V$ is taken by analysis functions. The interpolation problem then is the dual problem of synthesis where signal in $V$ is retrieved from the measurements in $\ell2$.

Sampling

The traditional Shannon’s sampling procedure is usually done as follows:

\[x(t) \rightarrow \text{Ideal Low Pass Filter} \rightarrow \text{Sample at regular interval} \rightarrow \{y_{k}\}\]

Alternatively, to see this mathematically, we have that:

\begin{align} y_k &= (x * \varphi) (k) , k\in \mathbb{Z} \\
&= \int_{\mathbb{R}} x(s)\varphi(k-s)ds \\
&= \int_{\mathbb{R}} x(s)\tilde{\varphi}(s-k)ds \\
&= \langle x(\cdot), \tilde{\varphi}(\cdot-k) \rangle_{L2} \\
&= \langle x, \tilde{\varphi}_{k} \rangle \end{align}

where \(\tilde{\varphi}_k(t) = \varphi(k-t)\)

is the measurement function of the sampling procedure.

Therefore, the sampling process can be looked-upon as

projection on shifted-verions of the time-reversal of LPF’s impulse response.

Synthesis (Interpolation)

The synthesis procedure (interpolation) is achieved by a weighted sum of the synthesis functions $\psi_k$ by the measurements $y_{k}$.

\[\hat{x}(t) = \sum_{k \in \mathbb{Z}} y_{k} \psi_{k} = \sum_{k \in \mathbb{Z}} y_{k} \psi(t-k)\]

Analysis+Synthesis can be simultaneously achieved using $sinc$

In general we have that \begin{equation} \hat{x}(t) = \sum_{k \in \mathbb{Z}} y_{k} \psi_{k} = \sum_{k \in \mathbb{Z}} \langle x(\cdot), \tilde{\varphi}(\cdot-k) \rangle \psi(t-k) = \sum_{k \in \mathbb{Z}} \langle x(\cdot), \varphi(k-\cdot) \rangle \psi(t-k) \end{equation}

In the case of where $\varphi(t) = sinc(t)$, which is a symmetic real function, both analysis and synthesis can be achieved by the same functional form as shown below:

\begin{equation} \hat{x}(t) = \sum_{k \in \mathbb{Z}} y_{k} \psi_{k} = \sum_{k \in \mathbb{Z}} \langle x(\cdot), sinc(\cdot-k) \rangle sinc(t-k) \end{equation}

The reason why we can use $sinc$ for both analysis and synthesis is shown below.

Interlude: Getting $sinc$ functions in Numpy

The numpy implementation of $sinc$ function has the following form:

\begin{equation} sinc(\omega t) = \frac{sin(\omega \pi t)}{\pi t} \end{equation}

The Fourier Transform is given by:

\begin{equation} \hat{f}(\xi) = \int_{-\infty}^{\infty} f(x) e^{-2\pi i x \xi} dx \end{equation}

As such, the DFT of $sinc$ function is given as

\begin{equation} \mathcal{F}[sinc(\omega t)] = \mathbb{1}_{[-\omega/2,\omega/2]} \end{equation}

Since our function are defined as arrays (they are discer-time in nature), the relevant fourier transform to invoke is the DFT. For numpy.fft, the DFT formulation for a length $n$ signal $x[k], k = 0,\ldots,n-1$ is given as

\[A_{k} = \sum_{m=0}^{n-1} a_m \exp\bigg\{-2\pi i \frac{mk}{n}\bigg\} \qquad k= 0,\ldots,n-1\]

Where the IDFT is given by: \(a_{m} = \frac{1}{n}\sum_{k=0}^{n-1} A_k \exp\bigg\{2\pi i \frac{mk}{n}\bigg\} \qquad m= 0,\ldots,n-1\)

Note: We are following the convention that the DFT is not normalized where the IDFT is normalized by $1/n$

Also Note: Most-likely, you will see ringing artifact of the fourier transform of $sinc$ and it is necesarrily caused by the sharp cutoff of the ideal low pass filter ($\mathcal{F}[sinc]$). To mitigate this problem, a hamming window can be applied prior to taking FFT.

dt = 1e-4 # sampling time step
Fs = 1./dt # sampling frequency
t = np.arange(-0.5,0.5,dt) # time vector
N = len(t) # number of sample-points

omega = 10.
x = np.sinc(omega*t)
w = np.hamming(len(x))
x2 = x*w
#x = np.sin(200.0 * 2.0*np.pi*t)
n = int(2**(np.floor(np.log2(N))+3)) # sample precision = nextpower2
sp = np.fft.fft(x,n=n)
sp2 = np.fft.fft(x2,n=n)
freq = np.fft.fftfreq(n,d=dt)

plt.figure(figsize=(15,6))
plt.subplot(121)
plt.plot(t,x,label=r'$sinc(t)$')
plt.plot(t,x2,label=r'$sinc(t)\odot W_{hamming}$')
plt.title(r'$sinc(\omega t), \omega=%.1f \ [rad^{-1}]$'%(omega),fontsize=20)
plt.xlabel(r'$t$',fontsize=15)

plt.subplot(222)
plt.plot(freq, abs(sp),label=r'$\mathcal{F}[sinc(t)]$')
plt.plot(freq, abs(sp2),label=r'$\mathcal{F}[sinc(t)\odot W_{hamming}]$')
plt.xlim([-20,20])
plt.title(r'$\|\mathcal{F}[sinc](\omega)\|$',fontsize=20)

plt.subplot(224)
plt.plot(freq, sp.real,label=r'$Re\{\mathcal{F}(sinc)\}$')
plt.plot(freq, sp.imag,'--',label=r'$Im\{\mathcal{F}(sinc)\}$')
plt.xlim([-20,20])
plt.legend()
plt.xlabel(r'$\omega [Hz]$',fontsize=15)
plt.tight_layout()
plt.show()

/Users/tingkailiu/anaconda3/envs/site/lib/python3.6/site-packages/mkl_fft/_numpy_fft.py:158: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  output = mkl_fft.fft(a, n, axis)

png

$\varphi_k, \psi_k$ :: ${sinc(\cdot-k)}$ form orthonormal basis

The reason why shifted $sinc$ functions can be used for both analysis and synthesis is clear if we consider the fact that they form an orthonormal basis for the bandlimited space $V$. We illustrate this fact by the following correlation plot between $sinc(t)$ and $sinc(t-k), k \in \mathbb{Z}$. Note that if the signals $\phi_{k}(t) = sinc(t-k)$ are indeed orthonormal, their innerproduct $<\phi_k,\phi_l>$ should be $\delta_{k,l}$.

delay_v = np.arange(-10,10,1) 
x_cor = np.zeros_like(delay_v)
for idx,delay in enumerate(delay_v):
    x_cor[idx] = np.correlate(np.sinc(omega*t),np.sinc(omega*(t-delay)))

    
x_cor = x_cor/np.sum(np.sinc(omega*t)**2) #normalize by energy of sinc
plt.figure(figsize=(15,5))
plt.stem(delay_v,x_cor)
plt.title(r'$\varphi = sinc$ form orthonormal basis')
plt.xlabel(r'Time Shift $k$')
plt.ylabel(r'$Corr(\varphi(t),\varphi(t-k))$')
plt.show()

png

Demo: Sampling & Interpolation with $sinc$

To put Shannon’s Sampling to test, we illustrate the sampling procedure on a randomly generated trigonometric polynomials \begin{align} x(t) = \frac{1}{\sqrt{T_{M}}} \sum_{n=-M}^{M}c_n \exp\bigg(j n \frac{\Omega}{M} x\bigg)
T_{M} = 2\pi M/\Omega \end{align}

The sampling procedure with sampling interval $T$ is \begin{align} y_k = (h \odot x)(t)\rvert_{t = nT} \end{align} while the interpolation is given by \begin{align} \hat{x} = \sum_{k} y_k \cdot sinc\big(\frac{t-kT}{T}\big) \end{align}

Create Input with gaussian random coefficients \begin{align} c_k = a_k + i\cdot b_k
a_k,b_k \sim \mathcal{N}(0,1) \end{align}

dt = 1e-3
Fs = 1./dt # sampling frequency
M = 10 # degree
Omega = 2*np.pi*15. # bandwidth
T_M = 2*np.pi*M/Omega # Period

t = np.arange(0.,T_M*5,dt)
#t = np.arange(0.,5.,dt)

a = np.random.rand(M)
b = np.random.rand(M)
c = a+1j*b

x = np.zeros((len(t),),dtype=np.complex128)
for n in range(M):
    if n == 0:
        x +=  c[n]
    else:
        x += c[n]*np.exp(1j*n*Omega/M*t)
        x += np.conj(c[n])*np.exp(1j*-n*Omega/M*t)
x = np.real(x)

N = len(t) # number of sample-points
n = int(2**(np.floor(np.log2(N))+1)) # sample precision = nextpower2
sp = np.fft.fft(x,n=n)
freq = np.fft.fftfreq(n,d=dt)

/Users/tingkailiu/anaconda3/envs/site/lib/python3.6/site-packages/mkl_fft/_numpy_fft.py:158: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  output = mkl_fft.fft(a, n, axis)

plt.figure(figsize=(15,5))
plt.subplot(121)
plt.plot(t,x)
plt.title(r'$x(t)$')
plt.xlabel(r't')
plt.subplot(122)
plt.plot(freq,abs(sp))
plt.xlim([-Omega/2/np.pi,Omega/2/np.pi])
plt.title(r'$\mathcal{F}[x]$')
plt.xlabel(r'$\omega [Hz]$')
plt.show()

png

Pre-filter with $sinc$

$sinc$ form an ideal low pass filter, therefore we can just take the FFT and truncate at desired frequency.

import scipy.signal
import copy

plt.figure(figsize=(15,8))
w_sinc = 10. # bandwidth of sinc in Hz

# LPF with sinc
sp = np.fft.fft(x)
freq = np.fft.fftfreq(len(sp),d=dt)
sp_hx = np.zeros_like(sp)
sp_hx[np.logical_and(freq>=-w_sinc,freq<= w_sinc)] = sp[np.logical_and(freq>=-w_sinc,freq<= w_sinc)]
hx = np.real(np.fft.ifft(sp_hx))
                         
ax1= plt.subplot(221)
ax2= plt.subplot(222)
ax3= plt.subplot(223)
ax4= plt.subplot(224)

ax1.plot(t,x)
ax3.plot(freq,abs(sp))
ax3.set_xlim([-25,25])
ax2.plot(t,hx)
ax4.plot(freq,abs(sp_hx),label=r'$\mathcal{F}[h \odot x]$')
ax4.legend(loc='upper right')
ax4.set_xlim([-25,25])

ax1.set_title(r'$x(t)$',fontsize=20)
ax1.set_xlabel(r'$t$',fontsize=18)
ax2.set_title(r'$(h * x)(t)$',fontsize=20)
ax2.set_xlabel(r'$t$',fontsize=18)
ax3.set_title(r'$\mathcal{F}[x]$',fontsize=20)
ax3.set_xlabel(r'$f [Hz]$',fontsize=18)
ax4.set_title(r'$\mathcal{F}[h * x]$',fontsize=20)
ax4.set_xlabel(r'$f [Hz]$',fontsize=18)
plt.tight_layout()
plt.show()

png

Sampling at Nyquist Rate

Since the pre-filtered signal has bandwidth w_sinc, we need to sampling at twice that frequency. Alternative, we set the sampling interval to be 1/(2*w_sinc)

T_sample = 1./(2.*w_sinc)
sample_interval = int(T_sample//dt)
y = hx[::sample_interval]  # sample amplitude
t_y = t[::sample_interval] # sample locations

plt.figure(figsize=(15,5))
plt.plot(t,hx,'k--',linewidth=1,label=r'$(h\odot x)(t)$')
plt.stem(t_y,y,markerfmt='C2o',label=r'$(h\odot x)(k), T = {0} [s]$'.format(T_sample))
plt.xlabel(r'$t$',fontsize=15)
plt.legend(loc='upper right',fontsize=15)
plt.title(r'Prefilter+Sampled @ $T = {0:.1e}$'.format(T_sample),fontsize=20)
plt.show()

png

Interpolation with $\varphi_k = sinc(t-k)$

x_rec = np.zeros_like(t)
d = t_y[1]-t_y[0]
for ak,tk in zip(y,t_y):
    x_rec += ak*np.sinc((t-tk)/(sample_interval*dt))


plt.figure(figsize=(15,5))
plt.plot(t[::sample_interval],y,'C2o',label=r'$(h\odot x)(k), T = {0} [s]$'.format(T_sample))
plt.plot(t,x_rec,label=r'$\hat{x}(t)$')
plt.plot(t,hx,label=r'$(h\odot x)(t)$')
plt.xlabel(r'$t$',fontsize=15)
plt.legend(loc='upper right',fontsize=15)
plt.title(r'Interpolation at Nyquist Rate with $\varphi_{k} = sinc(\omega (t-k))$',fontsize=20)
plt.show()

png

Comparison of Recovery at different Sampling Rate

The indices corresponding to Nyquist rate sampling interval is 50, we look at sampling intervals above and below this value.

sample_intervals = np.array([30,50,100,200])
colors = matplotlib.cm.Accent(np.linspace(0,1,len(sample_intervals)))

x_rec = np.zeros((len(sample_intervals),len(t)))
for i,T in enumerate(sample_intervals):
    y = hx[::T]  # sample amplitude
    t_y = t[::T] # sample locations

    for ak,tk in zip(y,t_y):
        x_rec[i,:] += ak*np.sinc((t-tk)/(T*dt))


plt.figure(figsize=(15,8))
ax1=plt.subplot(211)
ax2=plt.subplot(212)

for i,T in enumerate(sample_intervals):
    ax1.plot(t,x_rec[i,:],label=r'$\hat{x}(t)$ @ T=%.1e' % (T*dt),color=colors[i,:])
    ax2.plot(t,10*np.log10(abs(hx-x_rec[i,:])),color=colors[i,:])

ax1.plot(t,hx,label=r'$(h\odot x)(t)$',linestyle='--',linewidth=3)

ax1.set_xlim([min(t),max(t)])
ax2.set_xlim([min(t),max(t)])
ax2.set_xlabel(r'$t$',fontsize=15)
ax2.set_ylabel(r'RMSE [dB]',fontsize=15)
ax2.set_ylim([-100,50])
ax1.legend(loc='upper right',fontsize=15)
ax1.set_title(r'Interpolation Accuracy vs. Sampling Rates with $\varphi_{k} = sinc(\omega (t-k))$',fontsize=20)
plt.show()

/Users/tingkailiu/anaconda3/envs/site/lib/python3.6/site-packages/ipykernel_launcher.py:19: RuntimeWarning: divide by zero encountered in log10

png