Recall that a dirty image is formed using the following fourier relation between the spatial coherence and observed sky: \begin{equation} V_{measured}(u,v,w) = \left< G_p(\lambda,t)\left(\int_{sources}J_p(l,m,n,\lambda,t) I(l,m,n)e^{\frac{2\pi i}{\lambda}(ul+vm+w(n-1))}J_q^H(l,m,n,\lambda,t)\frac{dldm}{n}\right)G_q^H(\lambda,t) \right> \end{equation}
Provided that $||\Delta{w}||2\pi\left(\sqrt{1-\Delta{l}^2-\Delta{m}^2} - 1\right) << 1$ where $\Delta{x} = x_{max} - x_{min}$, we can take the inverse fast fourier transform (planar approximation) to get out a dirty image. However, we first have to resample the image (and if we don't sample at nyquest this has to include an anti-aliasing filter):
\begin{equation} V_{gridded}(u,v) = (V_{measured}\cdot S) * (C\cdot W)_{oversampled} \cdot (III\cdot II) \end{equation}Now we may take the inverse fast fourier transform to get the dirty image out: \begin{equation} V_{gridded}(u,v) \leftrightharpoons^\mathcal{\over{F}} I_{dirty}(l,m)\\ \end{equation}
Here the windowing function ($W$) is normally just a box ($II$). However: the ultimate goal is to limit any aliasing energy outside the box. Therefore we want to drive down the sidelobes of the fourier response using different windowing functions.
In [ ]:
import numpy as np
from matplotlib import pyplot as plt
from scipy import special as sp
from scipy import signal as sg
%matplotlib inline
In [ ]:
half_support = 20
full_support = 2 * half_support + 1
oversample = 63
beta = (7.43 - 2.39) / (5.0-2.0) * full_support - 0.366833
x = ((np.arange(0, full_support*oversample) - full_support*oversample/2))/float(oversample)
In [ ]:
def energy_outside_image(fft_filter,step=0.00001):
intg_lbound = (full_support/2 + 1)*oversample
return np.trapz(np.abs(fft_filter[intg_lbound:]),dx=step)*2 #symetric bound
In [ ]:
def energy_over_filter(fft_filter,step=0.00001):
return np.trapz(np.abs(fft_filter),dx=step)*2 #symetric bound
In [ ]:
sinc = np.sinc(x)
sinc /= np.sum(sinc)
fft_sinc = np.real(np.fft.ifftshift(np.fft.fft(np.fft.fftshift(sinc))))
plt.figure(figsize=(16, 6))
plt.title("BOX")
plt.plot(x,np.abs(sinc))
plt.xlabel("u")
plt.show()
plt.figure(figsize=(16, 6))
plt.title("BOX RESPONSE")
plt.xlabel("Image width")
plt.plot(x,np.abs(fft_sinc))
intg_value = x[(full_support/2 + 1)*oversample]
plt.plot([intg_value,intg_value],[np.min(fft_sinc),np.max(fft_sinc)],'r-')
plt.plot([-intg_value,-intg_value],[np.min(fft_sinc),np.max(fft_sinc)],'r-')
plt.show()
In [ ]:
sinc_hamming = np.sinc(x) * np.hamming(full_support*oversample)
sinc_hamming /= np.sum(sinc_hamming)
fft_sinc_hamming = np.real(np.fft.ifftshift(np.fft.fft(np.fft.fftshift(sinc_hamming))))
plt.figure(figsize=(16, 6))
plt.title("HAMMING")
plt.plot(x,np.abs(sinc_hamming))
plt.xlabel("u")
plt.show()
plt.figure(figsize=(16, 6))
plt.title("HAMMING RESPONSE")
plt.xlabel("Image width")
plt.plot(x,np.abs(fft_sinc_hamming))
intg_value = x[(full_support/2 + 1)*oversample]
plt.plot([intg_value,intg_value],[np.min(fft_sinc_hamming),np.max(fft_sinc_hamming)],'r-')
plt.plot([-intg_value,-intg_value],[np.min(fft_sinc_hamming),np.max(fft_sinc_hamming)],'r-')
plt.show()
In [ ]:
sqrt_inner = 1 - (2*x/float(full_support))**2
sinc = np.sinc(x)
kb = 1.0/full_support * sp.i0(beta * np.sqrt(sqrt_inner)) * sinc
kb /= np.sum(kb)
fft_kb = np.real(np.fft.ifftshift(np.fft.fft(np.fft.fftshift(kb))))
plt.figure(figsize=(16, 6))
plt.title("KAISER BESSEL")
plt.plot(x,np.abs(kb))
plt.xlabel("u")
plt.show()
plt.figure(figsize=(16, 6))
plt.title("KAISER BESSEL RESPONSE")
plt.plot(x,np.abs(fft_kb))
plt.xlabel("Image width")
intg_value = x[(full_support/2 + 1)*oversample]
plt.plot([intg_value,intg_value],[np.min(fft_kb),np.max(fft_kb)],'r-')
plt.plot([-intg_value,-intg_value],[np.min(fft_kb),np.max(fft_kb)],'r-')
plt.show()
In [ ]:
dx = 0.000001
aSinc=energy_outside_image(fft_sinc,dx)/energy_over_filter(fft_sinc,dx) * 100.0
aKBSinc=energy_outside_image(fft_kb,dx)/energy_over_filter(fft_kb,dx) * 100.0
aHammingSinc=energy_outside_image(fft_sinc_hamming,dx)/energy_over_filter(fft_sinc_hamming,dx) * 100.0
print "% Aliasing energy with rect window:",aSinc
print "% Aliasing energy with kb window:",aKBSinc
print "% Aliasing energy with hamming window:",aHammingSinc
In [ ]: