Generalized Sampling
In this notebook, we extend Shannon’s Sampling theorem to general Frame Theory settings.
Shift-Invariant-Space
We first extend the discussion in Bandlimited-Space to Shift-invariant-space:
\[V(\varphi) = \bigg\{s(x) = \sum_{k\in \mathbb{Z}} c_k \varphi(x-k): c\in \ell_2 \bigg\}\]Intuitively, $V(\varphi)$ is invariant over shift/time-delays:
\begin{equation} \forall f \in V(\varphi), f(x-k) \in V(\varphi), k \in \mathbb{Z} \end{equation}
Since shift in time-domain is equivalent to multiplication in fourier domain by an exponential, the bandwidth of the shifted signal remains unchanged, therefore
$L_{2} \cap B_{[-\Omega/2,\Omega/2]}$ is also a shift-invariant space
With this generalization, we can consider more general cases of bases functions, which we would discuss in the next sec tion.
Sampling and Riesz Basis
Consider analysis function $\varphi = sinc$ where $\tilde{\varphi}(x) = sinc(-x) = sinc(x) = \varphi(x)$:
\begin{align}
f(x) &= \sum_{k\in \mathbb{Z}} c_k \varphi_k \\
&= \sum_{k\in \mathbb{Z}} \langle f(\cdot), \tilde{\varphi}(\cdot-k) \rangle \varphi_k \\
&= \sum_{k\in \mathbb{Z}} \langle f(\cdot), \varphi(\cdot-k) \rangle \varphi_k \\
&= \sum_{k\in \mathbb{Z}} \langle f(\cdot), \varphi_k \rangle \varphi_k \\
\end{align}
It is clear that ${\varphi_{k} = sinc(t-k)}$ form an orthonormal basis for $L_{2} \cap B_{[-\Omega/2,\Omega/2]}$. The act of sampling then can be viewed as taking measurements of the original signal over this basis.
To invoke a similar measurement intuition for general shift-invariant-space, we introduce Riesz basis (see here for details) which can be linearly depdendent but still forms a stable and unambiguous representation of the original sigals. In another word, the basis functions \(\varphi_{n}\) need to satisfy the following condition:
The above definition of Riesz basis has the following implications:
- The projection onto the Riesz Basis is stable in that a small variation in the coefficients ${c_k}$ should lead to a small variation in the the function. The upper bound shown above keeps the $L^2$ norm of the function in check.
- Projection onto the Riesz Basis is isomorphic.
In addition to the above requirement, we would also want the model to
have the ability of approximating any input function as closely as desired by selecting a sampling step that is sufficiently small The above condition can be shown to be equivalent to the Partition of unity condition: \begin{align} \forall x \in \mathbb{R}, \sum_{k\in \mathbb{Z}} \varphi(x+k) = 1\\
y_{k} = \langle x(\cdot), \tilde{\varphi}(\cdot-k) \rangle, \tilde{\varphi}(t) = \varphi(-t) \end{align}
Reforumalating Norm-equivalence in Fourier Space
The first condition for Riesz Basis says that the norms are equivalent, suggesting that this projection from $L_{2}$ to $\ell_{2}$ onto ${\varphi_{k}}$ is isomorphic. This norm-equivalence condition, however, is hard to check because it requires checking $\forall c_k \in \ell_{2}$. An equivalent condition that is much easier to check is given as follows:
\[A \leq \sum_{k\in \mathbb{Z}} |\hat{\varphi}(\omega+2\pi k)|^{2} \leq B\]which says that the summation of energy of basis functions need to be bounded from both directions.
Note: This condition is equivalent to the one above because of the shift-invariant structure of the space we assumed.
Orthonormal Bases are Riesz Bases
As expected, all orthonormal bases are riesz bases. More importantly, it can be shown that a Riesz Basis is orthonormal if and only if:
\[A = B = 1\]We will now illustrate this fact with $sinc$ kernels. We use the following fact:
\[\varphi(t)\cdot \exp(j \omega_{0} t) \xrightarrow{\mathcal{F}} \hat{\varphi}(\omega+\omega_{0})\]Note that suppose that we change the bandwidth of the sinc kernel to be $sinc(\omega t)$ (i.e. $B_{[-\Omega/2, \Omega/2]}$), then the we can find it’s time-domain representation as
\[sinc(\omega t)\cdot \exp(j \omega_{0} \omega t) \xleftarrow{\mathcal{F}^{-1}} \hat{sinc}(\omega+2\pi k)\]Note: The $sinc$ function is passed through a hamming window before taking FFT to reducing the ringing artifacts.
Note: The amplitude of the fourier transform is not normalized and therefore is not 1. This is by convention.
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
w = np.hamming(len(t))
fft_n = int(2**(np.floor(np.log2(N))+1)) # sample precision = nextpower2
omega = 20.
freq = np.fft.fftfreq(fft_n,d=dt)
shift_range = np.arange(-3,4)
sp = np.zeros((len(shift_range),fft_n),dtype=np.complex128)
fig=plt.figure(figsize=(15,6))
plt.suptitle(r'$sinc$ satisfy Norm-Equivalence with A=B=1',fontsize=20)
ax1 = plt.subplot(121)
ax2 = plt.subplot(222)
ax3 = plt.subplot(224)
for k_i,k in enumerate(shift_range):
x = np.sinc(omega*t) * np.exp(-1j*2*np.pi*k*omega*t)
x *= w
sp[k_i,:] = np.fft.fft(x,n=fft_n)
ax2.plot(freq, abs(sp[k_i,:]),color=plt.cm.cool(k_i/len(shift_range)))
ax1.plot(t,x)
ax3.plot(freq,np.sum(abs(sp),axis=0))
ax1.set_title(r'Example Trace $sinc(t) \cdot \exp(j \omega_{0}t)$',fontsize=18)
ax1.set_xlabel(r'$t$',fontsize=15)
ax2.set_xlim([-100,100])
ax3.set_xlim([-100,100])
ax2.set_title(r'$\|\hat{sinc}(\omega+2\pi k)\|$',fontsize=18)
ax3.set_title(r'$\|\sum_{k} \hat{sinc}(\omega+2\pi k)\|$',fontsize=18)
ax3.set_xlabel(r'$\omega [Hz]$',fontsize=15)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
dt = 1e-4 # sampling time step
Fs = 1./dt # sampling frequency
t = np.arange(-10,10,dt) # time vector
N = len(t) # number of sample-points
w = np.hamming(len(t))
omega = 1.
shift_range = np.arange(-20,21)
x_all = np.zeros((len(shift_range),len(t)))
fig=plt.figure(figsize=(15,5))
ax1 = plt.subplot(111)
cs = plt.cm.Greys(np.linspace(0,1,len(shift_range)))
for k_i,k in enumerate(shift_range):
x_all[k_i,:] = np.sinc(omega*t-k)
ax1.plot(t,x_all[k_i,:],color=cs[k_i,:])
ax1.plot(t,np.sum(x_all,axis=0),'r',label=r'$\sum_{k}sinc(\omega t+k)=1$')
ax1.set_title(r'$sinc$ satisfy partition-of-unity',fontsize=20)
ax1.set_xlabel(r'$t$',fontsize=18)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
Other Riesz Basis Functions - B-spline
In practice, $sinc$ kernels are rarely used bacause it has inifinite temporal support. To mitigate this situation, other basis are used, including the Basis-Splines (B-Splines). The B-splines are pariticularly attactive for 2 reasons:
- N-th order B-spline ($\beta_{N}(t)$) can reproduce polynomials of maximum degree N
- Short (compact) support and excellent approximation properties
- Any function $\varphi(t)$ that reproduces polynomials of maxiumum degree N can be decomposed into a B-spline and a distibution $u(t): \int u(t) dt \neq 0$, i.e. $\varphi(t)=u(t) * \beta_{N}(t)$
The B-spline can be defined recursively as follows:
\(\beta_{0} = 1_{[-0.5,0.5]}\) \(\beta_{N}(t) = \beta_{0} * \beta_{N-1}(t), N\geq 1\) The temporal support for B-splines are $[-\frac{n+1}{2},\frac{n+1}{2}]$.
dt = 1e-4 # sampling time step
Fs = 1./dt # sampling frequency
t = np.arange(-3,3,dt) # time vector
N = len(t) # number of sample-points
w = np.hamming(len(t))
fft_n = int(2**(np.floor(np.log2(N))+1)) # sample precision = nextpower2
# setup FFT frequency
freq = np.fft.fftfreq(fft_n,d=dt)
# number of b-splines
bspline_n = 5
# color for plot
colors = matplotlib.cm.jet(np.linspace(0.,1.,bspline_n))
cs = plt.cm.cool(np.linspace(0.1,1,bspline_n))
# b-spline array, shape= (number, length)
beta = np.zeros((bspline_n,len(t)))
beta[0,np.logical_and(t>=-0.5, t<=0.5)] = 1.
# setup plots
plt.figure(figsize=(15,6))
ax1 = plt.subplot(121)
ax2 = plt.subplot(122)
ax1.plot(t,beta[0,:],color=cs[0,:],label= r"N=0")
ax2.plot(freq,abs(np.fft.fft(beta[0,:],n=fft_n)),color=cs[0,:],label= r"N=0")
for n in range(1,bspline_n):
beta[n,:] = dt*np.convolve(beta[0,:],beta[n-1,:],mode='same')
ax1.plot(t,beta[n,:],color=cs[n,:],label= r"N={0}".format(n))
ax2.plot(freq,abs(np.fft.fft(beta[n,:],n=fft_n)),color=cs[n,:],label= r"N={0}".format(n))
ax1.set_xlabel(r'$t$',fontsize=15)
ax1.set_title('First {0} Centered B-splines'.format(bspline_n),fontsize=20)
ax2.set_xlim([-5,5])
ax2.set_xlabel(r'$\omega [Hz]$',fontsize=15)
ax2.set_title(r'$\mathcal{F}[\beta_{N}]$',fontsize=20)
ax1.legend(fontsize=15)
ax2.legend(fontsize=15)
plt.show()
B-spline satisfy Riesz Condition
For this part, demonstrate that B-spline of degree 5 satisfy Riesz Condition, note that the support for $\beta_{5}$ is $[-3,3]$
import scipy.signal
beta_5 = scipy.signal.bspline(t, 5)
dt = 1e-4 # sampling time step
Fs = 1./dt # sampling frequency
t = np.arange(-10,10,dt) # time vector
N = len(t) # number of sample-points
w = np.hamming(len(t))
fft_n = int(2**(np.floor(np.log2(N))+1)) # sample precision = nextpower2
freq = np.fft.fftfreq(fft_n,d=dt)
shift_range = np.arange(-6,7)
sp = np.zeros((len(shift_range),fft_n),dtype=np.complex128)
x_all = np.zeros((len(shift_range),len(t)),dtype=np.complex128)
fig=plt.figure(figsize=(15,6))
plt.suptitle(r'$\beta_{5}$ satisfy Norm-Equivalence',fontsize=20)
cs = plt.cm.cool(np.linspace(0,1,len(shift_range)))
ax1 = plt.subplot(121)
ax2 = plt.subplot(222)
ax3 = plt.subplot(224)
for k_i,k in enumerate(shift_range):
x_all[k_i,np.logical_and(t>=k-3,t<k+3)] = beta_5
x_all[k_i,:] *= np.exp(-1j*2*np.pi*k*t)
sp[k_i,:] = np.fft.fft(x_all[k_i,:],n=fft_n)
ax2.plot(freq, abs(sp[k_i,:]),color=cs[k_i,:])
ax1.plot(t,x_all[k_i,:],color=cs[k_i,:])
ax3.plot(freq,np.sum(abs(sp),axis=0))
ax1.set_title(r'Shifted $\beta_{N}(t+k)$',fontsize=18)
ax1.set_xlabel(r'$t$',fontsize=15)
ax2.set_xlim([-5,5])
ax3.set_xlim([-5,5])
ax2.set_title(r'$\|\hat{sinc}(\omega+2\pi k)\|$',fontsize=18)
ax3.set_title(r'$\|\sum_{k} \hat{sinc}(\omega+2\pi k)\|$',fontsize=18)
ax2.set_xlabel(r'$\omega [Hz]$',fontsize=15)
ax3.set_xlabel(r'$\omega [Hz]$',fontsize=15)
plt.tight_layout(rect=[0, 0.03,1.,0.95])
plt.show()
B-spline satisfy Partition of Unity Condition
dt = 1e-4 # sampling time step
Fs = 1./dt # sampling frequency
t = np.arange(-10,10,dt) # time vector
N = len(t) # number of sample-points
w = np.hamming(len(t))
omega = 1.
shift_range = np.arange(-6,7)
x_all = np.zeros((len(shift_range),len(t)))
fig=plt.figure(figsize=(15,6))
ax1 = plt.subplot(111)
cs = plt.cm.Greys(np.linspace(0.1,1,len(shift_range)))
for k_i,k in enumerate(shift_range):
x_all[k_i,np.logical_and(t>=k-3,t<k+3)] = beta_5
ax1.plot(t,x_all[k_i,:],color=cs[k_i,:])
ax1.plot(t,np.sum(x_all,axis=0),'r')
ax1.set_title(r'$\beta_{5}(t)$ satisfy partition-of-unity $\sum_{k}\beta_{5}(t+k)=1$',fontsize=20)
ax1.set_xlabel(r'$t$',fontsize=15)
plt.show()
