Radio Astronomical Image Deconvolution Using Prolate Spheroidal Wave Functions

0.34 MB·English

by  Sarod Yatawatta

RADIO ASTRONOMICAL IMAGE DECONVOLUTION USING PROLATE SPHEROIDAL
WAVE FUNCTIONS
Sarod Yatawatta
Kapteyn Astronomical Institute, University of Groningen, Groningen,
and ASTRON, Dwingeloo,
1
The Netherlands.
1
0
2
n ABSTRACT Due to the noise floor, any extended source has finite support [4],
a In order to produce high dynamic range images in radio in- andtheregionofinterest(ROI)orthesupportofagivenextended
J terferometry, bright extended sources need to be removed with source might not be optimal for a given arbitrary basis.
minimal error. However, this is not a trivial task because the In order to tackle the problem of finding the optimal basis for
4
Fourier plane is sampled only at a finite number of points. The a given extended source, independent of the observed data, we
1
ensuing deconvolution problem has been solved in many ways, selectprolatespheroidalwavefunctions(PSWF)[5],[6].Asimilar
mainlybyalgorithmsbasedonCLEAN.However,suchalgorithms problem has been solved in magnetic resonance imaging [7] and
] that use image pixels as basis functions have inherent limitations we extend that result to radio interferometry in this paper. Unlike
M
and by using an orthonormal basis that span the whole image, we dataderived basisfunctions, PSWFbasiscanbe precomputed and
I can overcome them [1]. The construction of such an orthonormal can be reused for observations at different epochs. Furthermore,
. basis involves fine tuning of many free parameters that define the unlikeshapelets,wesufferlessfromartifactsoutsidetheROIwith
h
basisfunctions. Theoptimal basisfor agiven problem (or agiven minimalnumberofbasisfunctionsused.Infact,PSWFarealready
p
extendedsource)isnotguaranteed.Inthispaper,wediscusstheuse beingusedinradioastronomicalimagingtoconstructaregulargrid
o- of generalized prolate spheroidal wave functions as a basis. Given of sampling points inthe Fourier plane [8]. Werefer thereader to
r thegeometry (or the region of interest) of an extended source and [9] for similar applications of PSWFin geoscience.
t the sampling points on the visibility plane, we can construct the Notation: We denote vectors in bold lowercase and matrices in
s optimal basis tomodel thesource. Not only does thisgives us the bolduppercase.Thematrixtranspose,Hermitian,pseudoinverseare
a
minimum number of basis functions required but alsothe artifacts denoted by (.)T,(.)H and (.)† respectively. Theidentitymatrixis
[
outside the region of interest are minimized. given by I.

1
IndexTerms— Radioastronomy, Radiointerferometry,Decon- II. MATHEMATICAL FOUNDATIONS
v
volution
We present the basics of interferometric imaging and the use
0
of PSWF in this section. For a complete overview of radio
3 I. INTRODUCTION
interferometry, the reader is referred to [8].
8 Due to the increase in computing power, high dynamic range
2 imaginginradiointerferometry(onlylimitedbycalibrationerrors) II-A. Interferometric imaging
. isachievableandisessentialtoproducenovelscientificresults.One
1
0 ofthemainobstaclestoreachahighdynamicrangeisthefactthat v m
onlyafinitenumberofsamplesoftheFourier(visibility)planedata
1
is observed. Moreover, in earth rotation synthesis, these sampling
1
points are not on a regular grid. If the observed field of view
:
v containsbrightextendedsources,completeremoval(deconvolution)
i ofsuchsourcestorevealthefaintbackgroundremainsachallenge.
X The commonly used algorithm for such problems is CLEAN [2].
N
r However,asshownin[1](andreferencestherein),theusageofaset y
a of image pixels (clean components) in CLEAN based algorithms N
has limitations. b
An alternative approach to using clean components is to use
an orthonormal basis to model bright extended structure in the
observation. Once such a model is obtained, it can be subtracted N
x
from the visibilities to reveal the residual (or the fainter back-
ground). The most obvious method of deriving such a basis is to
use the observed data itself [3]. However, this has the drawback u l
(a) (b)
thatthederivationofthebasiscanonlybedoneoncethecomplete
observation is available (problematic for real time imaging etc.).
On the other hand, we can adopt any arbitrary orthonormal basis Fig. 1. (a) Sampling points (total N ) in the Fourier (visibility)
a
to a given observation, completely independent of the data. For plane. (b) Imageof N by N pixels, withthesupport (ROI)area
x y
instance, in [1] Gauss-Hermite polynomials (shapelets) were used shaded. The shaded area has N pixels.
b
to produce high dynamic range images. The drawback in such an
approachistheselectionoffreeparameters(suchasthenumberof
basisfunctions, andthe scale) cannot bedetermined inan optimal WeconsiderthevisibilityplanetobecomposedofN sampling
a
fashion. It requires experience of the user as well as a trial and points as in Fig. 1 (a). The image has N (= N ×N ) pixels.
x y
error approach to fine tune such a basis for a given observation. However,theROIhasonlyN pixels,correspondingtotheshaded
b
area in Fig. 1 (b). The ROI is determined by the support of the to find the mode vector m (size M by 1).
source structure and the noise floor [4]. 4) Find equivalent basis Pb= THP in the Fourier plane and
Let f(u ,v ) be the sampled visibility at the p-th point in the subtract the model frome the observed data z,
p p
visibilitey plane. Let us also denote the intensity of the q-th pixel
ontheimageasf(l ,m ).Thesetwoquantitiesarerelatedbythe r=z−Pm (8)
van Cittert-Zernikeqtheoqrem [8] and can be approximated by the eb
to get the residual data vector r.
Fourier transform for images with small support as:
ThecompleteproofofthederivationofPSWFisgivenin[7].For
Na−1 completeness,wegiveasketchofproofasfollows.Thenumerator
f(lq,mq)= X f(up,vp)ej2π(lqup+mqvp), q∈[0,N −1] (1) of (4) can be written as
e
p=0
kITpk2 (=a) kITTHpk2 =trace(pHTI ITTHp) (9)
N−1 b b b b
f(up,vp)= Xf(lq,mq)e−j2π(lqup+mqvp), p∈[0,Na−1]. (=b) trace(ζeHR−HTI ITeTHR−1ζ) e
e b b
q=0
Wecan represent (1) in vectorized form as
Here, (a) is obtained by the substitution p = THp and (b) is
f = Tf, f =THf (2) obtained by using ζ =△ Rp, where R=△ (TTH)1/2.eIt is easy to
where ef =△ [f(l ,m ),...e,f(l ,m )]T, simplify the denominator oef (4) by substituting Ib =I in (9) as
0 0 N−1 N−1
f =△ [f(u ,v ),...,f(u ,v )]T. kpk2 =trace(ζHR−HTTHR−1ζ)(=c)trace(ζHζ) (10)
0 0 Na−1 Na−1
e e e
ThematrixT (sizeN ×N) has on itsq-throw and p-thcolumn We used the fact that R−HTTHR−1 =I to obtain (c) in (10).
a
e−j2π(lpuq+mpvq).Notethatgivenf,anestimateforthetrueimage Using (9) and (10), we can rewrite (4) as
can be formed as f = THf but wee recover the true image only ζ = argmax trace(ζHR−HTI ITTHR−1ζ) (11)
if THT =I. In orbder to saetisfy this condition we need a regular b b
gridof sampling pointsinthevisibilityplane. Thisisnot thecase ζ, kζk=1
in reality for radio interferometry because the sampling points are Thesolutionto(11)isthelargesteigenvalue,eigenvectorpairofthe
determined by location of individual receivers and earth rotation. matrix K = R−HTI ITTHR−1. The dimension of K (N by
b b a
N ) maekes the computation of eigendecomposition proehibitively
II-B. Prolate spheroidal basis a
expensive. In order to reduce this, we apply the transform η =△
Our objective is to find a basis function p(l,m) (and its coun-
ITTHR−1ζ to Kζ =λζ to get
terpartinthevisibilityplanep(u,v))thatmaximizestheenergyin b
e
the ROI. The criterion for seleection can be written as IHTHR−2TI η=λη (12)
b b
|p(l,m)|2
λ= P(l,m)∈ROI . (3) which gives us the kernel K in (5) which is of dimension Nb by
|p(l,m)|2 N .
P(l,m) b
Inotherwords,λin(3)istheratiobetweentheenergyconcentrated II-C. Computational cost reduction
intheROIandthetotalenergyintheimage.Thehigher thevalue Under the assumption N ≫ N > N , the cost of computing
a b
of λ, we have lower sidelobes and artifacts outside the ROI. (5) involves finding the pseudoinverse (TTH)† (size N by N )
a a
Let the vectorized versions of p(l,m) and p(u,v), evaluated at whose rank is at most N. However, the rank of K will be at
thepixelsandvisibilitypoints,bepandp,resepectively. Then,we most N . Instead of using T, wedownsample the visibilitypoints
b
can rewrite (3) as e to construct a matrix T of size N by N (N ≥ N). This
λ= kITbpk2 (4) downsampling can be cbombined withdthe imagingd weights used.
kpk2 Therefore,(5)and(6)areevaluatedusingTinsteadofT.However,
when we calculate the residual in (8), webuse the full matrix T.
The selection of pixels that belong to the ROI from the full
image vector p is done by premultiplying by IT. Thus, IT is an III. EXAMPLE
b b
N ×N matrix, constructed by removing rows (corresponding to
b We take a LOFAR (http://www.lofar.org) test observation of
pixels outside the ROI) from an N ×N identity matrix.
Cygnus A, at a frequency of 213 MHz as an example. The image
Westatethemainresulthere,whichisadirectextensionof[7].
of the source Cygnus A is given in Fig. 2(a). The peak flux (not
The steps required for the removal of an extended source from a
normalized) is about 20 Jy and the total is about 10 kJy. The
given observation are:
objectiveistosubtractthissourcefromtheobserveddatatoseethe
1) Construct the kernel K (size Nb by Nb) as faint background sources. The ROI of this source, which is above
the noise floor is given in Fig. 2 (b). The image has N = 7552
K=ITbTH(TTH)†TIb (5) pixelsofdimension118by64andtheROIhasNb =1826pixels.
The observation lasted for about 8 hours and the original
2) FindtheeigendecompositionofKandselecttheeigenmodes
sampling points are shown in Fig. 3 (a). The number of sampling
with largest eigenvalues. If the i-th eigenvalue, eigenvector
points is N ≈15×106. We downsample the sampling coverage
pair of K is (λ ,η ), the i-th basis vector is a
i i (select only a subset with uniform probability) to get N =8931
d
1 whichisjustaboveN.Abetterapproachfordownsamplingwould
p = TH(TTH)†TI η (6)
i λ b i betocombinethiswithimagingweights(i.e.,selectingfewershort
i
baselines and more long baselines).
3) Represent theimageasthevector b(sizeN by1).Decom- Thedominant eigenvalues ofthekernel K(5)(dimension 1826
pose image b, using M basis vectors P=△[p ,...,p ] by 1826) is shown in Fig. 4. We see that approximately the first
0 M−1
100 eigenvalues are close to 1 while any eigenvalue beyond the
Pm=b, m=P†b (7) 200-th eigenmode is almost zero.
b
In Fig. 5, we have shown the PSWF basis vectors p corre-
sponding tothefirstfew largesteigenvalues. Notethatthesupport
of the basis vectors are entirely within the ROI, thus creating
minimal artifacts outside the ROI. We select the first 100 basis
functions M = 100 to construct the matrix P in (7). Using this,
we decompose the image in Fig. 2 (a) to get the mode vector m.
Using the equivalent basis inthe Fourier plane, to calculate the
residualasin(8).Oncetheresidualisobtainedwemaketheresidual
images as shown in Fig. 6. For comparison we have also shown
the residual image obtained using a shapelet basis function based
deconvolution.Bothmethodsgivearesidualnoiselevelofabout11
mJyfarawayfromCygnusA.Incontrast,withtraditionalCLEAN
based deconvolution, we get a residual noise of about 13 mJy.
Theshapelet based model usedabout 300basisfunctions ofmany
(a) (b)
scales. On the other hand, the PSWF model used only about 100
modes, which is much less.
Fig. 3. The sampling points in the Fourier plane (a) Original
coverage,foran8hourobservationwithabout15millionsampling
points. (b) Reduced coverage by downsampling by 1500 with
8931 sampling points. The downsampling is done with uniform
probability but it is also possible to combine this with imaging
weights, for instance by selecting fewer short baselines and more
long baselines.
1
0.9
0.8
(a) 0.7
0.6
λ0.5
0.4
0.3
0.2
0.1
0
0 50 100 150 200
Eigenvalue Number
(b)
Fig. 4. Eigenspectrum of the kernel K with the first 200 largest
Fig. 2. Cygnus A (a) Observed image at 213 MHz, deconvolved eigenvalues plotted.
using CLEAN. The peak flux is about 20 Jy and to total flux is
about 10 kJy. (b) The shaded area correspond to the ROI of 1826
pixels.
and Multichannel Signal Processing Workshop (SAM), Israel,
pp. 69–72, 2010.
IV. CONCLUSIONS [2] J.Hogbom,“Aperturesynthesiswithanonregulardistribution
ofinterferometerbaselines,”A&ASuppl.,vol.15,pp.417–426,
Wehavepresented theuseof prolatespheroidal wavefunctions
1974.
in radio interferometric image deconvolution and have demon-
[3] R. Levanda and A. Leshem, “Adaptive selective sidelobe can-
strated its feasibility by application to a real observation. We get
celler beamformer with applications in radio astronomy,” in
results comparable with existing techniques, but with fewer basis
proc. IEEE 26-th Convention of Electrical and Electronics
functionsand withlessartifactsoutsidetheROI.Futurework will
Engineers (IEEEI),Israel, 2010.
focus on widefield imaging and reducing the computational cost.
[4] D.Slepian,“Onbandwidth,”Proc.oftheIEEE,vol.64,no.3,
pp. 292–300, 1976.
V. ACKNOWLEDGMENTS
[5] D. Slepian and H. O. Pollak, “Prolate spheroidal wave func-
WethankWimBrouwandGerdeBruynforvaluablecomments tions, Fourier analysis, and uncertainty-I,” Bell Syst. Tech. J.,
and advise. vol. 40, no. 2, pp. 43–61, 1961.
[6] H. J. Landau and H. O. Pollak, “Prolate spheroidal wave
VI. REFERENCES
functions,Fourieranalysis,anduncertainty-II,”BellSyst.Tech.
[1] S. Yatawatta, “Fundamental limitations of pixel based image J., vol. 40, no. 2, pp. 65–84, 1961.
deconvolutioninradioastronomy,”inproc.IEEESensorArray [7] M. A. Lindquist, C. Zhang, G. Glover, L. Shepp, and
Q. X. Yang, “A generalization of the two dimensional prolate
spheroidal wave function method for nonrectiliner MRI data
acquisition methods,” IEEE Trans. on Image Proc., vol. 15,
no. 9, pp. 2792–2804, 2006.
[8] W. N. Brouw, “Aperture synthesis,” in Methods in Computa-
tional Physics, vol. 14, pp. 131–175, 1975.
[9] F. J. Simons and D. V. Wang, “Spatiospectral conventration
in the Cartesian plane,” Appl. Comput. Harmon. Anal. under
review, 2010.
Fig. 5. Some PSWF basis vectors corresponding to the largest
eigenvalues. The ROI is indicated by the dark curve.
(a)
(b)
Fig.6.ResidualsaftersubtractingCygnusAusing(a)prolatebasis
(b) shapelet basis. Note that in (b) there ismore unsubtracted flux
towards the center of the source while in (a) the residuals are
compact (point like). Both images have peak values of 5.5 Jy.