Fundamental Imaging Limits of Radio Telescope Arrays

VERSION 11 JUNE 2008 1 Fundamental Imaging Limits of Radio T elescope Arrays Stefan J. W ijnholds Studen t member , IEEE and Alle-Jan v an der V een F ellow , IEEE Abstract —The fidelity of radio astr onomical image s is generally assessed by practical e xperience, i.e. using rules o f thumb, although some aspects and cases hav e been treated rigorously . In this paper we present a m athematical framewo rk capable of describing the fund amental limits of radio astronomica l imaging problems. Alt hough the data model assumes a single snapshot observ ation, i.e. v ariations in time and frequency are not considered, this framework is sufficiently general to allow extension to synthesis observ ations. Using tools from statistical signal process ing and linear algebra, we discu ss the tractability of the imaging and deconv olu tion problem, the redistribution o f noise in the map by t he imaging and decon volution process, the cov ariance of the image values due t o propagation of calibration errors and thermal noise and the upper limit on th e number of sources tractable by self calibration. The combin ation of cov ariance of the image v alues and the number of tractable sources determines the effectiv e noise floor achi ev able in the imaging process. The effective noise provides a better figure of merit than d ynamic range sin ce it includes the spatial variations of the noise. Our results provide han dles f or i mpro ving the imaging perform ance by design of the array . Index T erms —radio astronomy , imaging, deconv olution, noise, dynamic range I . I N T RO D U C T I O N The radio astro nomical commun ity is curr ently building or developing a numb er o f ne w instruments such as the lo w frequen cy arra y (L OF AR) [1], the squar e kilo meter ar ray (SKA) [ 2] and the Mileura wid e field array (MW A) [3]. Imaging an d s elf calibration of th ese ra dio telescopes will be com putationally demandin g tasks due to th e large numb er of array elemen ts. Muc h research is th erefore fo cused on finding clev er sho rt-cuts to reduce the amount of processing required , such as w -projec tion [ 4] o r facet imagin g [5] and different variants of CLEAN [6]. The validity a nd quality of these m ethods is generally assessed by p ractical experien ce. Attempts to do a rigor ous analysis are do ne fo r som e aspects and cases [7]–[10], but m ost of the time rules of thumb are used. This paper presen ts th e first co mprehen si ve mathematical framework capable of describing the fun damental limits of radio as tronom ical imag ing prob lems. Th e data model used in this paper applies to sna pshot observations, i.e. variations in time a nd frequency are not consider ed. Howev er , using a multi-measur ement data mo del such as tho se in [11]–[1 3] This work was supported by the Netherla nds Foundatio n for Research in Astronomy (ASTR ON) and by NWO-STW under the VICI programme (DTC.5893). S.J. W ijnhol ds is with AST R ON, Dwingelo o, The Netherlan ds. A.-J. van der V een is with the Delft Unive rsity of T echn ology , Delft, The Netherlands. Email: wijnholds@a stron.nl, a.j.vande rvee n@tudel ft.nl it is straightforward to extend the data mod el to synthesis observation and still apply the fr amew ork de scribed herein . The resolution of the final image (or map ) is normally determined by the size and config uration of the array and th e spatial taper fun ction. Und er th e assumption that the sky is mainly em pty , i.e. that the image co ntains only a few sources, maps with high er resolutio n than pr edicted b y the array configur ation (superresolu tion) can be made using CLEAN. The Maximum Entropy Meth od (MEM) [14] impo ses a si milar constraint by aiming for a solu tion that is as feature less as possible. In the array pro cessing literature, superresolution is achieved by high resolution direction of arrival (DO A) estimation techniques such as MUSIC [15] and weighted subspace fitting [16 ], [1 7]. In all these ap proache s the go al is to disentangle the spatial re sponse of the array and the source struc ture, a proce ss called deconv olution. In section III we formulate imaging as an estimation problem, an appro ach called m odel based imaging , and o btain a n an alytic expression for its least squ ares solu tion that allows us to for mulate the deco n volution pr oblem as a m atrix inversion p roblem. This provides a po werful tool to assess the trac tability of the deconv olution problem and to dem onstrate the impact on the arr ay con figuration on the deco n volution p roblem an d the redistribution o f noise in the imag ing an d deconv olution process. The dynamic range of an im age is generally defin ed as the power r atio between the stron gest and th e weakest mean- ingful featur es in the map. In practice, the limitatio ns of an instrument are m ore conv eniently d escribed by the achievable noise floor in an im aging observation since the dynamic range strongly dep ends on the stren gth of the strongest source within the field -of-view an d b ecause the noise varies over the map. This noise floor is a combination of calibration errors, thermal noise a nd c onfusion noise. In th is p aper the ter m “effective noise” refers to th e net r esult of th ese constituents in the image p lane. In section I V an alytical expression s are de riv ed that describe the com ponents of the effecti ve no ise in terms of the covariance of the im age values, a con cept wh ich we will ref er to as image c ov ariance. The co nsequenc es of these expressions are illustrated with a f ew examples in section V. These examp les su ggest that the contr ibution of propagated calibration errors to the image cov ariance is considerab ly smaller than the contribution of thermal noise e ven if the calibration is d one on data with similar SNR. T hey also indicate that self calibration causes hig her cov ariance b etween source power estimates than pure im aging does. Notation : Overbar ( · ) denotes com plex conjug ation. The transpose operator is denoted by T , the complex c onjugate VERSION 11 JUNE 2008 2 (Hermitian) tran spose by H and the Moore-Pen rose p seudo- in verse by † . The expectation o perator is deno ted by E {·} , ⊙ is th e elemen t-wise matrix m ultiplication ( Hadamard produ ct), ( · ) ⊙ n is u sed to denote the elemen t-wise matrix expo nent with exponent n , ⊗ denotes the Kron ecker produ ct and ◦ is used to denote th e Khatri-Rao or colu mn-wise Kronec ker prod uct o f two m atrices. diag( · ) converts a vector to a diagonal matrix with th e elements of the vector placed on the main diag onal, vec( · ) converts a matrix to a vecto r by stacking the co lumns of the matrix and v ecdiag ( · ) converts the m ain diagon al of its argumen t to a co lumn vector . circulant( · ) creates a square circulant matrix by circularly shif ting th e entries of its vector argument to form its co lumns. ⊛ will be used to den ote circular conv olution of two vector s, i.e., for vectors of length n , ( x ⊛ y ) j = P n − 1 i =0 x i y j − i mo d n . For m atrices and vectors of compatible dimensions, we will frequen tly use the following p roperties: vec ( ABC ) = ( C T ⊗ A )vec( B ) (1) vec ( A diag( b ) C ) = ( C T ◦ A ) b (2) ( A ◦ B ) H ( C ◦ D ) = A H C ⊙ B H D (3) ( A ⊗ B )( C ◦ D ) = A C ◦ BD (4) I I . D A T A M O D E L Consider a phased array co nsisting o f p sensors ( an- tennas). D enote the baseb and outpu t signal o f the i th ar- ray elemen t as x i ( t ) an d d efine the array signal vector x ( t ) = [ x 1 ( t ) , x 2 ( t ) , · · · , x p ( t )] T . W e assume the presence of q source signals s k ( t ) impin ging on the array . The se are assumed to be m utually indepen dent i.i.d. Gaussian signa ls, and are stacked in a q × 1 vector s ( t ) . Like wise the sensor noise sig nals n i ( t ) ar e assumed to be mu tually ind ependen t Gaussian signals an d are stacked in a p × 1 vector n ( t ) . W e assume that the nar rowband condition ho lds [18]. W e can then describe, for the k th source signal, the phase delay differences over the p re cei ving elem ents du e to the propag ation geometry by a p -dimension al spatial signature vector a k . T he q spatial signature vectors are assumed to be kno wn (kno wn source locations and ar ray geometr y). The sensors are as sumed to ha ve the same direction de- penden t gain behavior which is described b y gain factors g 0 k tow ards the q source sign als re ceiv ed by the array . These ca n be collected in a m atrix G 0 = dia g ([ g 01 , g 02 , · · · , g 0 q ]) . The direction indep endent g ains and phases can b e described as γ = [ γ 1 , γ 2 , · · · , γ p ] T and φ = [e j φ 1 , e j φ 2 , · · · , e j φ p ] T respec- ti vely , with cor respond ing d iagonal matrix fo rms Γ = diag( γ ) and Φ = diag ( φ ) . W ith these definitions, the ar ray signal vector can be describ ed as x ( t ) = ΓΦ q X k =1 a k g 0 k s k ( t ) ! + n ( t ) = GAG 0 s ( t ) + n ( t ) (5) where A = [ a 1 , · · · , a k ] (size p × q ) and G = ΓΦ . The signal is sampled with period T and N sam - ple vectors are stacked into a data matrix X = [ x ( T ) , x (2 T ) , · · · , x ( N T )] . The covariance matrix o f x ( t ) is R = E { x ( t ) x H ( t ) } and is estimated by b R = N − 1 XX H . The number o f samples N in a snapsho t observation is equal to the p roduct of bandwidth and in tegration time an d typically ranges from 10 3 (1 s, 1 kHz) to 1 0 6 (10 s, 10 0 kHz) in radio as tronom ical applicatio ns. Likewise, the source signa l covariance Σ s = diag( σ s ) wher e σ s = [ σ 2 s 1 , σ 2 s 2 , · · · , σ 2 sq ] T and the noise covariance m atrix is Σ n = diag( σ n ) where σ n = [ σ 2 n 1 , σ 2 n 2 , · · · , σ 2 np ] T . Th en the m odel for the covari- ance matrix fo r a snapsho t observation R based o n (5) is R = GA G 0 Σ s G H 0 A H G H + Σ n . (6) If the d irectional response o f the antenn as is kno wn, G 0 can be absor bed in A . If G 0 and Σ s are both unk nown, we can introd uce Σ = G 0 Σ s G 0 H = diag([ | g 01 | 2 σ s 1 , · · · , | g 0 q | 2 σ sq ]) = diag( σ ) (7) with real valued elemen ts σ = [ σ 2 1 , σ 2 2 , · · · , σ 2 q ] T . W e m ay then restate (6) as R = GAΣA H G H + Σ n . (8) The i th ele ment of th e sensor ar ray is located at r i = [ x i , y i , z i ] T . Th ese positions can be stacked in a matrix R = [ r 1 , r 2 , · · · , r p ] T (size p × 3 ). The position of the k th source can be d enoted by the unit vector l k = [ l k , m k , n k ] T . The source po sitions can be stacked in a matrix L = [ l 1 , l 2 , · · · , l q ] T (size q × 3 ). The spatial signatur e matr ix A can thus b e described by A = exp  − j 2 π λ RL T  (9) where th e exponen tial fun ction is ap plied elem ent-wise to its argument. In the re mainder of th is p aper we will spe cialize to a planar array having z i = 0 for conv enience of presentation but w ithout loss of g enerality . I I I . I M AG I N G A N D D E C O N VO L U T I O N A. Beam fo rming versus model based imaging The imaging proce ss transfor ms the cov ariances of the received signals (called visibilities in radio astrono my) to an image of the source structure with in the field-of- view of the receivers. I n array processing terms, it can b e d escribed a s follows [11]. T o determine the p ower of a signal rece i ved from a particu lar dir ection ( l , m, n ) , a w eight vector w = ( a † ) H = exp  − j 2 π λ R [ l, m, n ] T  † H = 1 p exp  − j 2 π λ R [ l, m, n ] T  (10) is assigned to the array signal vector x ( t ) . The oper ation y ( t ) = w H x ( t ) is ge nerally called beamform ing an d can be regarded as a spatially matched filter . Equation (10) represents the m ost basic b eamforme r that assumes the presenc e of only a single so urce and on ly c orrects the sign al delays due to the array geometry . These weigh ts can b e adapted to correct the complex g ain d ifferences b etween the receiving elements G derived from calibration measuremen ts [ 19], nullin g of interfering sou rces [12] and spatial tapering of the arr ay [20]. VERSION 11 JUNE 2008 3 m l −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 0.2 0.4 0.6 0.8 1 1.2 m l −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 1.2 Fig. 1. ( a ) Image obtained by normal imaging without decon vol ution as in (11), showing the sources and thei r side lobe patterns. ( b ) Image obtained by model base d imagi ng as in (22), which estimate s the power at ev ery pixel simultane ously , resulti ng in a decon volved image sho wing only the sources without the array response . The imag e v alue at ( l , m, n ) is equal to th e expected output power of the beamform er wh en poin ted into that direction, and can be co mputed directly from the array covariance matrix b R as b i ( l , m, n ) = w H b Rw . (11) For weights defined as in (10), this is k nown as dir ect F ourier transform imagin g . T o create an image, w is scanned over all relev ant ( l, m, n ) . The re quired weights can be stacked into a single matrix W . Since w H b Rw = ( w ⊗ w ) H vec ( b R ) , we can stack all ima ge values in a single vector b i and wr ite b i = ( W ◦ W ) H vec ( b R ) . (12) If we on ly want to image a t the source location s, we h av e W = 1 p A . A typical model assumption is that there is a source present at every pixel location , in which case b i B F = 1 p 2 ( A ◦ A ) H vec( b R ) . (13) This is the classical dirty image . Let us assum e mo mentarily that G = I a nd Σ n = 0 . Inserting th e data mo del (6), or vec ( R ) = ( A ◦ A ) σ , into (13) gives i B F = E { b i B F } = 1 p 2 ( A ◦ A ) H ( A ◦ A ) σ = 1 p 2 ( A H A ⊙ A H A ) σ (14) This shows that the dirty image is not eq ual to the true sour ce structure. T o under stand the p hysical me aning of this term, consider th e produ ct a H i a j , where the indices i and j refer to the r espectiv e column s o f A . Using (9) th is can be written explicitly as a H i a j = exp  − j 2 π λ R l i  H exp  − j 2 π λ R l j  = p X n =1 exp  j 2 π λ r T n ( l i − l j )  . (15) The physical in terpretation of the inn er p roduct b etween the two spatial signatu re vector s is th at it m easures th e sensitivity of the ar ray to sign als com ing from directio n l j while the array is steered towards l i . The pro duct A H a j thus d escribes the array sensitivity for all d irections o f in terest stacked in L when pointe d to l j . It therefor e provide s the array voltage response or ar ray voltage beam p attern centered aro und l j , b V ( l j ) = A H a j . (16) W ith A defined as in ( 9), this shows that th e voltage beam pattern is just th e Fourier transfor m of the spatial weighting function resulting from the array co nfiguration and the weight- ing of th e ar ray elem ents. Th e corr esponding power beam pattern can be ca lculated as b P ( l j ) = b V ( l j ) ⊙ b V ( l j ) = A H a j ⊙ A H a j . (17) The factor A H A ⊙ A H A in (14) can thu s b e interpre ted as a conv o lution by the Fourier tran sform of the spatial distribution of baseline vectors, which is known as the array beam pattern or dirty bea m [21]. This e ffect is illustrated in Fig. 1( a ). This image is the result o f a simu lated observation with an 8 × 8 half wa velength spaced (i.e., spatially Nyquist sampled) 2D u niform rectangu - lar ar ray (URA). The grid of imag e values on th e sky is taken such that the first Nyqu ist zo ne is app ropriately samp led. T he underly ing source model conta ins four sou rces at grid poin ts ( − 0 . 33 , − 0 . 6 , 0 . 73) , ( − 0 . 2 , − 0 . 6 , 0 . 7 7) , (0 . 6 , − 0 . 2 , 0 . 77) and (0 . 87 , 0 . 2 , 0 . 46 ) respectively an d σ = [1 , 0 . 6 , 1 . 3 , 0 . 1 ] T . This source and array con figuration will be used th rough out this paper unless stated other wise. The map in Fig. 1( a ) clearly shows the se four (o r three, if o ne regard s the two sources on neighbo ring grid points as a sing le extended so urce) b eing conv o lved with the arr ay beam pattern. Follo wing a model based approach , the deco n volution prob- lem can be f ormulated as a maxim um likelihood (M L) esti- mation problem , that shou ld p rovide a statistically e fficient estimate of th e parameter s. Since a ll signals are assumed to be i. i.d. Gaussian sign als, the deriv ation is standard and the ML estimates are o btained by minim izing the negative log- likelihood function [22] b σ = argmin σ  ln | R ( σ ) | + tr  R − 1 ( σ ) b R  . (18) VERSION 11 JUNE 2008 4 It d oes no t seem possible to solve this minim ization prob - lem is closed fo rm, but a weig hted least squ ares covariance matching approach is k nown to lead to estimates that are, f or a large numb er of samp les, equivalent to ML estimates and therefor e asympto tically efficient [22]. The pro blem can thus be ref ormulated as b σ = argmin σ w w w W c  b R − Σ n  W c − W c GAΣA H G H W c w w w 2 F = argmin σ w w w  W c ⊗ W c  vec  b R − Σ n  −  W c GA  ◦ ( W c GA ) σ n w w w 2 F (19) The solution is given by b σ =  W c GA  ◦ ( W c GA )  †  W c ⊗ W c  vec  b R − Σ n  (20) Optimal weighting is provided by W c = R − 1 2 . Since radio as- tronom ical sources are generally very weak with th e stron gest source in th e field having an in stantaneous SNR in the orde r of 0.01 , we can introdu ce the appro ximation R ≈ σ 2 n I for an array o f identical elemen ts f or convenience of nota tion. This reduces (20) to b σ =  GA ◦ GA  † vec  b R − Σ n  . ( 21) One may argue that th is re quires one to kn ow where the sources are bef ore doing the imagin g. T his is gen erally so lved by simu ltaneously estimating the so urce location s an d sourc e powers. Although the C LEAN algor ithm h as not yet been fully analyzed, it can be r egarded as an iterativ e procedur e to do this [11]. It is in structiv e, howe ver , to use (21) for imaging by estimating the p ower on e very image p oint (pixel), i.e., by assuming a data mod el with a sou rce p resent at every p ixel. W e can simplify (21) by r eplacing th e Moore-Penr ose pseudo- in verse by the left pseudo- in verse, to o btain the image vector b i =   GA ◦ GA  H  GA ◦ GA   − 1 × ×  GA ◦ GA  H vec  b R − Σ n  =  A H Γ 2 A ⊙ A H Γ 2 A  − 1 × ×  GA ◦ GA  H vec  b R − Σ n  . (22) The first factor in this eq uation repr esents the d econv o lution operation . It is theref ore convenient to intr oduce th e deco n vo- lution matrix M = A H Γ 2 A ⊙ A H Γ 2 A =   A H Γ 2 A   ⊙ 2 . This provides a powerful check on the sam pling of the image plane. If the im age plane is oversampled , i.e., if too many im age points ar e defined , this matrix will be sing ular . This p roperty demonstra tes that h igh resolution imag ing is only possible if a limited numbe r o f sources is p resent, i.e. , if th e n umber of sources is mu ch sma ller than the numb er of resolution elements in the field -of-view . Th e cond ition nu mber of th e deconv olution matrix, which provides a measure on the mag- nification o f measurement noise, is discussed in more detail in Sec. III-C. This mo stly empty field-of -view is com monly assumed in astronom ical imaging and this assump tion is on e of the reasons why CLEAN and MEM work in pr actice. Fig. 1( b ) shows the image o btained by apply ing (22) to the 8 × 8 URA. Comparison with the image obtain ed using (1 3) clearly shows the effectiv eness o f the mod el based imaging appr oach in supp ressing th e array beam pa ttern. B. Noise redistrib u tion If imag ing is don e with out deconv olution b y using (13), the thermal noise adds a constant value to all imag e values. This can b e illustrated by assum ing that E n b R o = Σ n , i.e. by assuming that the im age is completely d ominated by thermal noise. The expected value of the imag e then become s i B F = 1 p 2  A ◦ A  H vec ( Σ n ) = 1 p 2  A ⊙ A  H σ n = 1 T σ n p 2 1 (23) where we used the fact that all elements of A ha ve u nit amplitude. Th is equation describes an image w here a ll values are equal to the average therm al no ise per baseline. If the ima ging pro cess in volves deconv olution, the result is described by (22). For simplicity we will assume that we have an ar ray of iden tical elem ents, so that we can set G = I . Furthe r , to illustrate the effect, we m omentarily omit th e correction by Σ n in (22). In this case, the expected value of the imag e is i =  A H A ⊙ A H A  − 1  A ◦ A  H Σ n =    A H A   ⊙ 2  − 1  A ⊙ A  H σ n =    A H A   ⊙ 2  − 1  1 T σ n  1 . (24) In this case, the hom ogeneity o f the the rmal noise d istribution in th e ma p d epends o n the row sums of    A H A   ⊙ 2  − 1 being constan t. If this is true, the mod el based imag e using (22) is analogou s to th e beamfor med image based on (13). A special case is the situation in which the co lumns of A are orthon ormal. Otherwise, the structur e is mor e co mplicated. Th is is il- lustrated in Fig. 2 which compar es the noise distribution in the im age p lane of the 8 × 8 URA b y assuming R = 0 . 1 I with the correspo nding image for a five armed array , ea ch arm bein g an e ight element half wa velen gth spaced Un iform Linear Array (ULA) . The impact of the r edistribution of noise can be redu ced by estimating the re cei ver no ise powers and subtracting these estimates fro m the array cov ariance matrix as described by (22). In most astrono mical imaging algor ithms, the autoco rrelations are g enerally ignored completely th us effecti vely introducin g a small negativ e system noise since the autocorr elations repr esent the power sum o f the source signals and the noise. VERSION 11 JUNE 2008 5 m l −1 −0.5 0 0.5 −1 −0.5 0 0.5 3 3.5 4 4.5 5 x 10 −4 m l −1 −0.5 0 0.5 −1 −0.5 0 0.5 3 3.5 4 4.5 5 x 10 −4 Fig. 2. ( a ) Im aging with decon voluti on using an 8 × 8 half wave length spaced array for a Nyquist sampled image assuming R = 0 . 1 I (an empty field with only the rmal noise). ( b ) Imaging result for a fi ve armed array , eac h arm being an eigh t element half wa vel ength spaced ULA. C. Deco n volution matrix co ndition number The deco n volution matrix M not only cau ses a redistribu- tion of noise over the map, but also determines whether the deconv olution is a well cond itioned problem. If the deconv o- lution matrix is no t inv ertible, the proble m is ill-po sed and additional con straints are req uired to ob tain a uniq ue solution. Different cho ices fo r these constrain ts or even the rigo r with which they are applied , lead to different imag ing results for CLEAN and MEM b ased on the same data. In some cases, this may ev en lead to different interpretation of the final maps [14]. These pr oblems arise d ue to over -interpre tation of the d ata by allo wing fo r more image points ( parameters) than can b e justified by the data. In these situation s, the cond ition num ber of the de con volution ma trix will be infinitely large. Even if the deconv olution m atrix is invertible, its condition numbe r may be unaccep tably high in view of the SNR of the data: the conditio n n umber is a measu re fo r the magn ification of measuremen t noise [ 23]. Th e con dition nu mber thus provid es a p owerful diagnostic too l to assess th e feasibility of th e deconv olution prob lem at hand. It is instru ctiv e to analyze a half wa velength s paced 1D ULA with ide ntical eleme nts, i.e. with G = I , sampling the sky on a regular g rid. In this case A represents a Fourier tran sform mapping th e spatial frequencies on th e sky to the spatial samples describing th e ele ctromagn etic field over the array aperture. As demon strated in th e previous section , these spatial frequen cies will be conv o lved in th e imag ing process with the Fourier tra nsform o f the arr ay ap erture tap er or voltage beam pattern, which can be easily calcu lated fo r l = 0 : b V ( 0 ) = A H a ( 0 ) = F T  1 p 0 n − p  . (25) Here F T d enotes the Fourier transfo rm, n is the total n umber of image points, p is the n umber of elem ents in the array and 0 p and 1 p denote p × 1 vectors c ontaining ze ros and one s respectively . The corresp onding power beam patter n is b P ( 0 ) = b V ( 0 ) ⊙ b V ( 0 ) = F T  1 p 0 n − p  ⊛  1 p 0 n − p  . (26) If th e co lumns o f A are ordered such that they describ e the array respon se vectors for the regula rly spaced DOAs starting with l 1 = 0 , it is easily seen that M = A H A ⊙ A H A = circulant ( b P ( 0 )) , (27) i.e., that th e deconv olution ma trix for a 1D ULA e quidistantly sampling th e image plane is a circulant matrix . Since M is a circulant matr ix, its eig en values λ = [ λ 1 , λ 2 , · · · , λ n ] T are giv en b y the Fourier transfor m o f b P ( 0 ) [24], or λ = F T ( b P ( 0 )) = F T  F T  1 p 0 n − p  ⊛  1 p 0 n − p  =  1 p 0 n − p  ⊛  1 p 0 n − p  (28) since F T ( · ) = F T − 1 ( · ) for real symmetric function s. For Hermitian matrices, the co ndition numb er κ is g i ven by the ratio o f the largest and smallest eige n value, i.e. κ = λ max /λ min [24]. If the ima ge plane is Ny quist samp led, n = 2 p − 1 and λ = [1 , 2 , · · · , p − 1 , p, p − 1 , · · · , 2 , 1] T . (29) In this case th e cond ition num ber o f M is κ = λ max λ min = p 1 = p, (30) thus M is inv ertible. The deco n volution problem is th erefore well-posed and has a uniqu e solutio n. If the image plane is unde rsampled with n < 2 p − 1 samples, then λ =  p − n − 1 2 , · · · , p − 1 , p, p − 1 , · · · , p − n − 1 2  T (31) and κ = 2 p 2 p − ( n − 1) . The deconv olution prob lem in itself is thus well-p osed and h as a un ique solution. Howe ver , from Fourier th eory we kn ow that aliasing effects may occu r du e to und ersampling . VERSION 11 JUNE 2008 6 0.05 0.1 0.15 0.2 0.25 0.3 0.35 10 0 10 5 10 10 10 15 10 20 ∆ l condition number 8 × 8 2D ULA five armed ULA Fig. 3. This plot shows the condition number of the decon volution m atrix as functio n of the image resolution for the 8 × 8 half wav elength spaced array and the fiv e armed array with each arm being an eight element half wa vel ength spaced ULA. If the image plane is oversampled with n > 2 p − 1 samples, then λ = [ · · · , 0 , 1 , 2 , · · · , p − 1 , p, p − 1 , · · · , 2 , 1 , 0 , · · · ] T (32) and κ = ∞ . In this c ase, the deconv olution problem is ill- posed an d thus no t solvable witho ut introdu cing addition al bound ary con ditions to constrain the pr oblem. This analy sis sho ws that, for a 1D ULA, the condition number slowly incr eases up to Nyqu ist sampling o f the image plane and then jumps to infinity . Since a URA is just the 2D an alog of a 1D ULA, this b ehavior is also expected for the 8 × 8 URA introduced ear lier . This co njecture is confirmed in Fig. 3 which shows the con dition numbe r of the deconv olution m atrix as fun ction of im age resolution . This figure a lso shows the co rrespond ing curve for the fiv e ar med array intro duced earlier to d emonstrate the impac t of less regular and sparser sampling o f the arr ay aperture. Althoug h the ar ray diameter is nearly twice as large, it does n ot provide twice the resolution du e to sparser samplin g of the ap erture plane. Th is plot also demonstrates that a less regular array also may have a less strict cu t-off: the transition of the con dition number from small values to infinity is a g radual one. For array processing proble ms, this m eans that the user shou ld decide which value of the condition n umber (or no ise enhan cement) is still accep table. Regularization is commo nly used to av oid u nin vertability of matrices. In radio astronom ical imag ing wh ere m ost sources have a low SNR, this would lead to imp erfect d econv o lution causing the weakest so urces in th e field to be dr owned in the imperfectly removed a rray respon se patter n of the strongest sources. Howe ver, se veral fo rms o f implicit regularization have been stud ied to handle special c ases like strong inter ference [25]. I V . E FF E C T I V E N O I S E Equation (22) shows that calibr ation and im aging are strongly cou pled. Knowledge o f the instru mental parameters is required to obtain th e p roper image. People have appro ached this prob lem in two ways. I n the first approach calibration and imaging are treated as sep arate steps, i. e., the instrumen tal parameters ar e estimated first b y a calibration measuremen t and consecutively ap plied to the actual m easurement data. The second approach is self calib ration wh ich regards the estimation of instrumen tal and image par ameters as a single parameter estimation pr oblem [13], [26]–[ 28 ]. In either case the ach ie vable d ynamic range is limited by the combinatio n o f estimation erro rs, thermal noise and con fusion noise. T ogeth er , they determine the effective noise in the image which need n ot be homog eneous over the field of interest. In this sectio n a n umber of analytica l expressions are derived that d escribe these contributions in terms of the d ata mo del presented in Section II. The im plications will be d iscussed in Section V. A. Noise in self calibrated images In self calibration th e instrumental and im age parameters a re estimated simultaneo usly . Self calibr ation based o n th e data model presented ab ove can thus be d escribed as simultan eous estimation of the om ni-direction al com plex gains, the app arent source p owers, th e sou rce location s and the receiver noise powers, i.e. , of a para meter vector θ = [ γ 1 , · · · , γ p , φ 2 , · · · , φ p , σ 2 2 , · · · , σ 2 q , σ 2 n 1 , · · · , σ 2 np , l 2 , · · · , l q , m 2 , · · · , m q ] T . In this parameter vector, φ 1 and σ 2 1 are omitted beca use they are set to constants for the problem to b e identifiable. Ind eed, the restriction σ 2 1 = 1 is imposed by th e fact that G and Σ share a commo n factor, while the first constrain t is required since o ne can o nly measure th e gain ph ases with r espect to some referen ce, h ere achieved by setting φ 1 = 0 . 1 Similarly to (1 9), the parameter s are o btained by solving b θ = ar gmin θ k W c  b R − GAΣA H G H − Σ n  W c k 2 F , (33) where G , A , Σ and Σ n are all f unctions of θ and W c = R − 1 2 ≈ 1 σ n I as argued earlier . The minimum v ariance f or an u nbiased estimator is given by the Cram ` er-Rao Bound (CRB). The CRB on th e erro r variance for any u nbiased estimator states that the covariance matrix C θ of the p arameter vector θ satisfies [30] C θ = E   b θ − θ   b θ − θ  T  ≥ 1 N J − 1 , (34) where J is the Fisher in formation matrix ( FIM). For Gaussian data models J can be expressed as (e .g. [ 31]) J = F H  R − 1 ⊗ R − 1  F (35) 1 In [29] it is shown that P p i =1 φ i = 0 is the optimal constraint for this problem. This constrain t has the disadv antage that the location of the phase referenc e is not well defined . Furthermore, the choice for the constra int used here simplifies our analysis in combinati on with the constraints required to uniquel y identify the source location s and the appa rent source powe rs. VERSION 11 JUNE 2008 7 where R is the data covariance matrix an d F is th e Jacobian ev alu ated at the tru e values of the parameter s, i.e., F = δ vec( R ) δ θ T    θ . (36) For th e self c alibration scenario, th e Jacobian can be parti- tioned into six parts fo llowing the structure of θ : F = [ F γ , F φ , F σ , F σ n , F l , F m ] . (37) By substitution of (8) in (36), it follows dire ctly that the first four com ponents ca n be expressed as F γ =  GR 0 Φ  ◦ I + I ◦  GR 0 Φ  (38) F φ = j  GR 0 G  ◦ I − I ◦  GR 0 G  I s (39) F σ =  GA  ◦ ( GA )  I s (40) F σ n = I ◦ I (41) where R 0 = AΣA H and I s is a selection ma trix o f ap pro- priate size equ al to th e iden tity matrix with its first co lumn removed so that the deriv ati ves with respect to φ 1 and σ 2 1 are omitted. If the re cei ver no ise p owers of all p elem ents are the sam e, the expression for F σ n giv en in (41) shou ld be replaced by F σ n = v ec( I ) (42) For the last two comp onents o f th e FIM, derivati ves of vec ( R ) with respect to the sou rce po sition coord inates are required . Let ( x i , y i ) be the coord inates of the i th array element, and in troduce G x = diag([ x 1 , x 2 , · · · x p ] T ) G (43) G y = diag([ y 1 , y 2 , · · · y p ] T ) G (44) then these co mponen ts c an be conveniently written as F l = − j 2 π λ  GA ◦ G x A − G x A ◦ GA  ΣI s (45) F m = − j 2 π λ  GA ◦ G y A − G y A ◦ GA  ΣI s (46) These equations show t hat the entries of the Jacobian related to deriv ati ves with respect to the l - and m -coordinate s of the sources are prop ortional to the x - an d y -coordinates of the array elements respectively . Th e physical interpr etation of this relatio n is that a plan e wa ve prop agating along the coordin ate ax is o f the coordinate to be estimated provides a more u seful test signal to estimate the sour ce location th an a signal pr opagating pe rpendicu lar to this axis. The precedin g equations allow us to compu te C θ . Th e variance of the estimated image values, i.e. th e no ise on the image values due to estimation inaccu racy , is g i ven by the diagona l of the sub-b lock C σσ of this matrix, following th e partitioning o f θ . In general C σσ is n ot a diagona l matrix. The other e ntries in th is sub-b lock describe the way in which the noise on the pixels are co rrelated am ong themselves—this is associated with false structu res. B. Pr opagation o f calibration err ors If the instrumental parameter s are extracted f rom separate calibration data, th e m inimum v ariance o n th ese estimated values is given by the CRB on th e instru mental pa rame- ters in the calib ration experimen t, C θ , where now θ = [ γ 1 , · · · , γ p , φ 2 , · · · , φ p , σ 2 n 1 , · · · , σ 2 np ] T . With th is cho ice for θ the results f or F γ , F φ and F σ n derived ea rlier in (38), (39) and (41) can be used assuming that the calibratio n measuremen t adheres to the same data model. The propa gation of the calib ration errors to the image is describ ed by cov ( i ) =  ∂ i ∂ θ T  C θ  ∂ i ∂ θ T  T . (47) W e thus need to derive ∂ i / ∂ γ T , ∂ i /∂ φ T and ∂ i /∂ σ T n . The derivati ve o f the imag e values to γ k is defin ed as ∂ i ∂ γ k = ∂ ∂ γ k M − 1  GA ◦ GA  H vec ( R − Σ n ) (48) where M and G d epend on γ . Applyin g the formula fo r th e deriv ati ve of an inverted m atrix with respect to one of its elements [24], this can be rewritten as ∂ i ∂ γ k = = − M − 1  ∂ ∂ γ k M  M − 1  GA ◦ GA  H + + M − 1 γ k  e − j φ k E kk A ◦ GA + GA ◦ e j φ k E kk A  H ! × vec ( R − Σ n ) (49) where E kk is the elementar y matrix with all its entries set to zer o excep t element E kk which is set to 1. Inserting the vectorized version of (8) in (4 9) and rem oving the Khatri-Rao produ cts, we o btain ∂ i ∂ γ k = − 2 γ k (2 − γ k ) M − 1 Re n A H k : A k : ⊙ A H Γ 2 A o σ (50) W e have introd uced the notation A k : = E kk A , i.e. A k : has only z ero valued entries excep t on the k th ro w where the elements are equal to the cor respondin g elements of A . The goal of th is de riv atio n is to obtain an expr ession for ∂ i /∂ γ T . W e will thus have to stack the expre ssion for ∂ i / ∂ γ k in a single matrix. Th is is facilitated by introduc ing a k as the k th row of A and rewriting (50) as ∂ i ∂ γ k = − 2 γ k (2 − γ k ) M − 1 Re  a T k 1 T ⊙ A H Γ 2 AΣa H k  (51) where 1 d enotes a vector of one s o f approp riate size. By stackin g all vector ∂ i /∂ γ k in a single matrix, we thus obtain ∂ i ∂ γ T = − 2 M − 1 Re  A T ⊙ A H Γ 2 AΣA H  (2 I − Γ ) Γ . (52) VERSION 11 JUNE 2008 8 The corr espondin g resu lt for φ k can b e derived in a similar way , so we on ly present the main step s. ∂ i ∂ φ k = = ∂ ∂ φ k  GA ◦ GA  † vec ( R ) = M − 1  ∂ ∂ φ k GA ◦ GA  H   GA ◦ GA  H σ  = M − 1  − je − j φ k Γ A k : ◦ GA + je j φ k GA ◦ ΓA k :  H ×   GA ◦ GA  H σ  . (53) Removal o f the Khatri- Rao products by redu cing them to Hadamard pro ducts g i ves ∂ i ∂ φ k = − 2 γ 2 k M − 1 Im n A H k : A k : ⊙ A H Γ 2 A o σ . (54) Note th at this term h as th e same form as the first term in (50), so it c an be rewritten in a similar way . T his gives ∂ i ∂ φ k = − 2 γ 2 k M − 1 Im  a T k 1 T ⊙ A H Γ 2 AΣa H k  (55) and theref ore ∂ i ∂ φ T = − 2 M − 1 Im n A H ⊙ A H Γ 2 AΣA H o Γ 2 (56) Finally , the partial de riv ative o f the image values with respect to ( σ 2 nk ) is given by ∂ i ∂ ( σ 2 nk ) = ∂ ∂ ( σ 2 nk ) M − 1  GA ◦ GA  H vec ( R − Σ n ) = M − 1  GA ◦ GA  H vec ( − E kk ) = − M − 1  | GA | ⊙ 2  H vecdiag ( E kk ) . (57) Therefo re ∂ i ∂ σ T n = − M − 1  | GA | ⊙ 2  H . (58) If Σ n = σ 2 n I this red uces further to ∂ i ∂ σ n = − M − 1  | GA | ⊙ 2  H 1 . (59) The partial deriv ati ves as well a s the CRB [19] con tain terms in volving A H A , o ften weighted by the gains o f the receiving eleme nts. Given the ph ysical inter pretation o f th is factor d iscussed in sectio n I II, this sugg ests th at the error patterns intr oduced in the im age by calibration err ors follow the structu res in th e dirty image. Th is is confirm ed by the example in section V. Since the CRB is in versely pro portion al with N , which is equ al to the prod uct of bandwidth and integration time, the imag e cov ariance due to calibration erro rs decreases prop ortional to band width and in tegration time . C. Therma l n oise In this section we derive an expression for the covariance of the imag e values due to the therma l n oise in the data. W e will there fore assume that p erfect knowledge of the ther mal noise p ower Σ n is available to a void confu sion between the thermal noise c ontribution an d th e con tribution of p ropaga ted estimation errors. T he covariance o f th e image values is b y definition given by cov ( i ) = = E   vec  b i  − vec ( i )   vec  b i  − vec ( i )  H  = E n  GA ◦ GA  †  vec  b R − Σ n  − vec ( R − Σ n )  ×  vec  b R − Σ n  − vec ( R − Σ n )  H ×  GA ◦ GA  † H o . (60) This shows that und er the assum ption that perfect knowledge on Σ n in b R is av ailable, Σ n drops out. Furthe rmore, ( GA ◦ GA ) † can b e moved outside the expectation op erator, since it contains no estimated values. There fore cov ( i ) = M − 1  GA ◦ GA  H cov ( R )  GA ◦ GA  M − 1 . (61) For Gau ssian data mod els cov ( R ) = 1 N  R ⊗ R  , (62) and we find that cov ( i ) = = 1 N M − 1  GA ◦ GA  H  R ⊗ R   GA ◦ GA  M − 1 . This can be rewritten using Kro necker and Khatr i-Rao prod uct relations as cov ( i ) = 1 N M − 1  GA ◦ GA  H  R GA ◦ R GA  M − 1 = 1 N M − 1   A H G H R GA   ⊙ 2 M − 1 . (63) Finally , substituting the data model pr esented in (8) we get cov ( i ) = = 1 N M − 1   A H Γ 2 AΣA H Γ 2 A + A H Γ 2 Σ n A   ⊙ 2 M − 1 (64) It is interesting to note that for an array having G = I a diagona lization of A H A does no t only ensure a ho mogen eous noise distribution over th e map after the deco n volution o per- ation as demonstrated in Sec. III-B, but also diag onalizes th e image cov ariance d ue to therma l n oise, thu s ensur ing that the noise on the pixels is unc orrelated. Th e Gram matrix A H A describes the amoun t o f linear indepen dence (o r orthogo nality) of the direction of arr i val vectors within the field o f view of the array , which can be visualized as th e arr ay b eam pattern. This observation th erefore suggests that an array with a low side lob e pattern does not only provid e good spatial separ ation between source signals, but also g i ves small covariance b etween image values after deconv olution. VERSION 11 JUNE 2008 9 D. Confusion noise The contributions to the ef fectiv e noise from calibr ation errors an d th ermal noise scale in versely w ith the n umber of samples N , which is eq ual to th e p roduct of ban dwidth and integration time. Th is implies that, th eoretically , these sources of image n oise ca n b e red uced to ar bitrarily low lev els. In practice , th e rad io astronom ical array will detect more sou rces with every redu ction of the noise in th e map. At some point, the nu mber of detected sour ces beco mes larger than the numbe r of re solution elements in the image, which will tur n the map into one blur of sources. The maximum density of discernable sources is the classical co nfusion limit and relates to the resolutio n of th e image. In terms of self calibr ation, to ha ve m ore detectable sources requires more source parameters to describe the source m odel. At some point, th e self calibration p roblem becomes ill-p osed. W e will refer to this as the self calibr ation con fusion limit. Although the exact limit depen ds on the minutiae of the array and source configuration, we can easily comp ute an up per limit on the tractab le n umber of sources based on the argumen t that the n umber of un knowns should be smaller than the numb er of equatio ns. The data mo del provides a relation between the parameters and the data. For a p element arr ay , the covariance matrix contains p 2 indepen dent real values, so the data mode l can be regard ed as p 2 indepen dent equations. Solving for the direction ind ependen t comp lex g ains requ ires 2 p − 1 real valued parame ters, estimation of the receiver element noise powers requir es another p par ameters and the q sou rces ar e described by 3 q − 1 param eters ( the ap parent sou rce powers relativ e to the first source and two co ordinates per sour ce). The self calibr ation prob lem is th erefore c onstrained by p 2 ≥ 2 p − 1 + p + 3 q − 1 = 3 p + 3 q − 2 (65) implying that q ≤ p 2 − 3 p + 2 3 . (66) The spatial Nyquist sampling with the 8 × 8 URA allows an im age g rid of 15 × 15 = 22 5 image values. This r esolution was confirm ed by th e cond ition nu mber an alysis presented in Fig. 3. Howe ver, the upper lim it based on the an alysis above f or a 64 e lement array is 1 302. Th e mismatch between this upper limit and the a ctual numb er of u niquely solvable image values can be attributed to the redun dancy in the ar ray . Due to this redund ancy the cross-corr elations of many an tenna pairs provide the same sp atial informa tion instead of providing additional inform ation on the spatial structure of the sky . In terms of the argumen t leading to Eq. (66), there is linear depend ence between th e equations and therefore the nu mber of eq uations th at can be used to solve parameters is red uced. The 5-arme d array perf orms much better in th is r egard. E .g., for p = 40 the up per limit on the num ber of sources given by (66) is 494 . Since √ 494 ≈ 2 2 , we can thu s form an image grid of 22 × 22 points, thus provid ing a resolutio n of ∆ l ≈ 0 . 09 . Figure 3 shows that the condition num ber f or this array goes to num erical infinity a t ∆ l ≈ 0 . 0 8 , showing that the 5 -armed array approach es its theoretical self calibr ation confusion l imit. This example illustrates that if processing po wer is cheap com- pared to an tenna hardware, a non- redund ant array should b e image index image index 50 100 150 200 50 100 150 200 −10 −9 −8 −7 −6 −5 −4 image index image index 50 100 150 200 50 100 150 200 −10 −9 −8 −7 −6 −5 −4 Fig. 4. ( a ) The logarit hm (base 10) of the image cov ariance matrix due to calib ration errors and ( b ) due to the measurement noise on the same color scale. The former is almost ev erywhere two orders of magnitude lower . preferr ed over a re dundan t array if the confusion limit should be reached without in troducin g an ill-p osed de con volution problem . In this section we addressed the classical co nfusion limit for p ure imag ing p roblems and the self c alibration confusion limit in self calib ration pr oblems. Th is type of confu sion is often called source confu sion as opposed to side l obe confusion which ref ers to blur ring of the imag e by side lobe leftovers introdu ced in the CLEAN process. In the analysis of this paper, side lobe confusio n is pa rt of the decon volution p roblem and is thus in trinsically in cluded in the analysis of calibration error and therm al no ise p ropagatio n, an d d oes not need to be addressed sep arately . Sou rce conf usion does require a separate treatment because it inv olves the source d ensity distribution as function of sou rce brightn ess. V . I M P L I C A T I O N S A. Thermal no ise vs. pr opagated ca libration err ors W e comp are the image covariance d ue to calibratio n error s to the imag e covariance d ue to the noise on the d ata in a simulation. The calibration paramete rs are calculated from a separate d ata set with the same data model and integratio n VERSION 11 JUNE 2008 10 20 40 60 80 100 120 10 −6 10 −5 10 −4 parameter index variance self calibration calibration observation Fig. 5. The CR B fo r th e omni -direct ional complex gain amplit udes (paramete rs 1 through 64) and phases (parameters 65 through 127) for a separat e calibration observa tion and the self calibra tion approach. time. W e comp uted th e CRB for the 8 × 8 URA using the relations presented in Section IV -A for the simultaneou s solution of th e omn i-directiona l complex gains, G , an d the system n oise power , σ 2 n , wh ich was assumed to be the same for all a rray elements, a ssuming a sh ort term integration over N = 16 384 samples. This CRB was u sed in (47) to compu te the imag e covariance m atrix du e to calibration err ors. The magnitud es o f the matrix entries are shown in Fig. 4( a ), in a log scale. The image covariance matrix due to the noise in the measuremen t was calculated using (64) a nd is shown in Fig. 4( b ). Th is shows that the covariance of th e im age values due to the calibr ation er rors is m ore conce ntrated at the source locations tha n the covariance du e to the system n oise, but is generally m ore than two order s of magn itude lower . T hese results indicate that the calibratio n er rors only represent a minor con tribution to the total effecti ve noise, even when the calibration m easurements are do ne on the same (short) time scale s u sing sources o f th e same streng th, i.e. when the calibration measuremen t is similar to the actu al measu rement. B. Calibration ob servations vs. self calibration In th e pr e vious section we discussed th e situation in which the ar ray is calibra ted in a separate measuremen t. This scheme requires an extremely stable instrumen t. In most pr actical applications, the calibra tion is therefore do ne on the sam e data that is also used to provide the final image (self calibra tion). It is interesting to see how th ese scen arios c ompare. T o this end, we used the relationship s p resented in Section IV -A to compute the CRB for simultaneou s estimation of the omni- directional co mplex g ains, the app arent source po wers, th e source locations and the system no ise power f or the 8 × 8 URA. Figure 5 co mpares the CRBs for th ese two cases. Th e expected covariance of the gains an d pha ses in th e self calibration experim ent is higher since mor e parameters have to be estimated simu ltaneously . The behavior of the CRB o n the T ABLE I C OVA R I A N C E O F S O U R C E P OW E R E S T I M A T E S . self cali bratio n inde x 2 3 4 2 0 . 266 × 10 − 4 0 . 287 × 10 − 4 0 . 022 × 10 − 4 3 0 . 287 × 10 − 4 1 . 237 × 10 − 4 0 . 048 × 10 − 4 4 0 . 022 × 10 − 4 0 . 048 × 10 − 4 0 . 007 × 10 − 4 separat e calibration inde x 2 3 4 2 0 . 280 × 10 − 4 0 . 072 × 10 − 4 0 . 006 × 10 − 4 3 0 . 072 × 10 − 4 0 . 772 × 10 − 4 0 . 012 × 10 − 4 4 0 . 006 × 10 − 4 0 . 012 × 10 − 4 0 . 005 × 10 − 4 phases in th e self calibrated observation (slop ed u pwards for increasing parameter ind ex) can be explaine d by th e interaction between the source param eters and the gain p hases comb ined with the choice of the phase referen ce elemen t in the co rner of the ar ray . T ab le I shows, fo r each o f the two cases, the c ov ariance matrices of the app arent source powers, i.e., the variance of the image values at the locations of the sou rces. The scaling factor ambig uity between G and Σ in the self calibration case is reso lved by putting σ 2 1 = 1 to constrain the problem, and theref ore only th e covariance values o f the oth er thre e sources is tabulated. For the case with a separate calibration observation th e covariance m atrix was extracted fro m the sum of th e image covariances due to calibratio n er rors and system noise. The results in the tab le indic ate that the variance o f the source power estimates in both cases are co mparable , although the source po wer estima tes ar e slightly be tter when g ain calibration data is av ailable f rom a sepa rate mea surement. The covariance values f ound for a sepa rate calibration stag e are much lower than the correspo nding values for self calibratio n. This sug gests that pur e imaging is mor e capable of separating source sig nals f rom d ifferent dir ections than self calibrated imaging. V I . C O N C L U S I O N S In this paper we presented an an alytic solution fo r snapshot imaging including deconv olution based on a d ata model (mea- surement e quation) for th e anten na signal covariance matrix or visibilities. The presented com prehensive framework is sufficiently flexible to enable extension o f this analysis to synthesis observations, since the data mo del f or a synthesis observation has the same form [11]–[1 3]. This framework allowed us to make th e first com plete rigorou s assessment o f the effecti ve n oise floor , which is the combined ef fect o f pr op- agated calib ration errors, th ermal noise and source co nfusion, in the image in ter ms of the covariance of the imag e values. Our simulations for a 2D unifo rm r ectangular array indicate that the effect of pro pagated calibration error s is strongly concentr ated at the source locations b ut is considerably s maller than the thermal noise at oth er ima ge p oints. Th e results also suggest that if the instrumen t is sufficiently stable, a separ ate calibration step is to be prefe rred over a self calib rated image since it allows better source separation in the imaging p rocess. VERSION 11 JUNE 2008 11 The effects of decon v olution can be describ ed by a decon v o- lution matrix th at descr ibes the amount of linear indepe ndence (ortho gonality) of the spatial signature vectors weighted by the actual gains of th e receiving elements. A d iagonal d econv o lu- tion matrix not o nly ensures the best possible spatial separ ation between th e sour ces, but also ensures a hom ogeneo us no ise distribution over the map. This po ses the question wh ether this matrix can be diago nalized by array design or by applying approp riate weig hts to the array elements. Sin ce this factor is related to the ar ray b eam pattern, the latter is equ iv alen t to finding we ights that supp ress th e side lobe patter ns at least in the d irection o f other sources, which sug gest that techn iques like Robust Capon Beamforming should provide the requested weighting [ 25]. The cond ition n umber of the d econv o lution matrix can be u sed to assess the quality of the solu tion to the deconv olution prob lem. Compared to a redu ndant array (ULA, URA), a n array without redund ant element spacings provides mu ch be tter possibilities to approach the max imum n umber of so lvable image po ints for a fixed n umber of a ntenna elements, thereb y allowing th e system to reach the theoretical self calibratio n confusion limit. R E F E R E N C E S [1] J. D. Bregman, “LOF AR Approaching the Critica l Design Revie w , ” in Pro ceedin gs of the XXV IIIth General Assembly of the International Union of Radio Science (URSI GA ) , New Delhi , India, Oct . 23-29 2005. [2] P . Hall, “The Square Kilometer Array: an Internati onal Engineering Perspecti ve, ” Experimental Astro nomy , v ol. 17, no. 1, pp. 5–16, 2004. [3] C. J. Lonsdale, R. J. Cappello, J. E. Salah , J. N. He witt, M. F . Morales, L. J. Greenhill , R. W ebste r and D. Barn es, “The Mileura W idefield Ar- ray , ” in Proce edings of the XXVIIIth Gener al Assembly of the Interna tio nal Union of Radio Sci ence (URSI GA) , New Del hi, India, Oct. 23-29, 2005. [4] T . J. Cornwell, K. Golap and S. Bhatnaga r , “ w -projection: A New Algorith m for Wi de Field Imaging with Radio Synthe sis Arrays, ” in ADASS XIV , ser . Astronomical Society of the Pacific Confer ence Series, vol. 347, 2005. [5] R. A. Perle y , “W ide Field Im aging II: Imaging with Non-coplan ar Arrays, ” in Synthesis Imaging in Radio Astrono my , ser . Astronomical Society of the Pacific Conference Series, Richard A. Perley , Fred- eric R. Schw ab, Alan H. Bridle , Ed. , vol. 6, 1994, pp. 139–165 . [6] T . Cornwell, R. Braun and D. S. Briggs, “Decon volution, ” in Synthesis Imagin g in Radio Astr onomy II , ser . Astronomica l Society of the Paci fic Conferen ce Series, G. B. T aylor , C. L. Carilli and R. A. Perle y, Ed., vol. 180, 1999, pp. 151–170. [7] U. Schwarz, “Mathemati cal-sta tistical Description of the Iterati ve Beam Removi ng T echnique (Method CLE AN), ” Astrono my & Astrop hyics , vol. 65, pp. 345–356, 1978. [8] S. T an, “An Analy sis of the Propert ies of CLEAN and Smoothness Stabili zed CLE AN—Some W arnings, ” Monthly Notices of the Royal Astr onomica l Society , vol. 220, pp. 971–1001, 1986. [9] S. R. Kulka rni, “Self-Noise in Interferomete rs: Radio and Infrared, ” Astr onomica l Journal , v ol. 98, no. 3, pp. 1112–1130, Sept. 1989. [10] S. J. W ijnholds, “Self-noi se in full sky LOF AR images, ” in Nederlandse Astr onomen Confer entie (NA C) , Ameland , The Nethe rlands, 10-12 May 2006. [11] A. Leshem and A. van der V een, “Radio astronomical ima ging in the presence of s trong radio interfere nce, ” IEE E T r . Information Th. , vol. 46, no. 5, pp. 1730–1747, Aug. 2000. [Onlin e]. A va ilabl e: ftp:// cas.et.tude lft.nl/pub/allejan- docs/it99.ps.gz [12] A.-J. va n der V een, A. Leshem and A.-J. Boonstra, “ Array Signal Processing for Radio Astronomy , ” Experi mental Astr onomy , v ol. 17, no. 1-3, pp. 231–249, 2004 . [13] S. v an der T ol, B. J ef fs, and A. van der V een, “Self calib ration for the L OF AR radio astronomical array , ” IEEE T r . Signal P r ocessin g , vol. 55, no. 9, pp. 4497 –4510, Sept. 2007. [Online]. A vail able: http:/ /ens.e w i.tudel ft.nl/pubs/jeffs06t sp.pdf [14] R. Narayan and R. Nityananda , “Maximum Entrop y Image Restora tion in Astronomy, ” A nnual R ev iew of Astr onomy & As tr ophysics , no. 24, pp. 127–170 , 1986. [15] R. O. Schmidt, “Multiple Em itter L ocati on and Signal Paramete r Es- timatio n, ” IEEE T rams. Antennas and Prop agat ion , vol. AP-34, no. 3, Mar . 1986. [16] M. V iber g and B. Ottersten, “Sensor Array Processing Based on Sub- space Fitting, ” IEEE T rans. Signal Pr ocessing , vol. 39, no. 5, pp. 1110– 1121, May 1991. [17] M. V iber g, B. Otterste n and T . Kailath, “Detecti on and Estimation in Sensor Arrays Using W eigh ted Subsp ace Fitting, ” IE EE T rans. Signal Pr ocessing , v ol. 39, no. 11, pp. 2436–2448, Nov . 1991. [18] M. Zatman, “How narrow is narro wband, ” IEE Pr oc. Radar , Sona r and Navig . , vol. 145, no. 2, pp. 85–91, Apr . 1998. [19] S. J. Wijn holds and A. -J. Boonstra, “A Multisource Calibrat ion Method for Phased Array Radio T ele scopes, ” in 4 th IEEE workshop on Sensor Array and Multichan nel Proc essing (SAM) , W altham (MA), USA, 12- 14July 2006. [20] S. J. W ijnho lds, “Redu cing the impact of station le ve l spatial filtering limitat ions, ” in SKA Calibr ation & Imaging W orkshop (calim) , Cape T own, South Africa, 4-6 Dec. 2006. [21] A. Thompson, J. Moran, and G. Swenson, Interfer ometry and Synthesis in Radio Astr onomy , 2nd ed. John W iley & Sons, Inc., 2001. [22] B. Ottersten, P . Stoica and R. Roy, “Cov arian ce Matching E stimation T echnique s for Array Signal Processing Applicat ions, ” Digital Signal Pr ocessing , A Revi ew J ournal , vol . 8, pp. 185–210, July 1998. [23] G. Golub and C. va n Loan, Matrix Computations . B altimore , MD: Johns Hopkins Uni versity Press, 1984. [24] T . K. Moon an d W . C. Sti rling, Mathemati cal Methods and Algorit hms for Signal Proce ssing . Upper Saddl e Riv er , New Jerse y: Prentice Hall, 2000. [25] S. va n der T ol and A. J. van der V een, “Applicat ion of Robust Capon Beamforming to Radio Astronomica l Imaging, ” in IEEE Int. Conf. on Acoustics, Speech and Signal Pr oc. (ICASSP) , Philadelphi a (P A), Mar . 2005. [26] T . Cornwell an d E. B. Fomalont, “Sel f-Cali bration , ” in Synthe sis Imaging in Radio Astronomy , ser . Astronomical Society of the Pa cific Conference Series, R. A. Perley , F . .R. Schwab and A. H. Bridl e, Ed. BookCrafte rs Inc., 1994, v ol. 6. [27] B. P . Flanagan and K. L. Bell, “Array Self-Cal ibrati on with Large Se nsor Position E rrors, ” in IEEE Internatioal Conf er ence on Acoustics, Speech and Signal Pr ocessing (ICASSP) , 1999. [28] T . J. Pearson and A. C. S. Readhead, “Image Formatio n by Self- Calibra tion in Radio As tronomy , ” A nn. Rev . Astron. Astr ophys. , v ol. 22, pp. 97–130, 1984. [29] S. J. W ijnhol ds and A. J. v an der V een, “Ef fects of Parame tric Con- straints on the CRLB in Gain and P hase Estimation Problems, ” IEEE Signal Pr ocessin g Letters , vo l. 13, no. 10, pp. 620–623, Oct. 2006. [30] S. Kay, Fundament als o f Statistical S ignal Pr ocessing: Estimation Theory . Engle wood Cliffs, Ne w Jersey: Prent ice Hall, 1993, vol. 1. [31] P . Stoica, B. Ottersten, M. V iberg, and R. Moses, “Maximum Likeli hood Array Processing for Stochast ic Coherent Sources, ” IEEE T ransactions on Signal Pr ocessi ng , vol. 44, no. 1, pp. 96–105, Jan. 1996. PLA CE PHO TO HERE Stefan J. Wijnholds (S’2006) was born in The Netherl ands in 1978. He recei ved M. Sc. degrees in Astronomy and Applied Physics (both cum laude) from the Univ ersity of Groningen in 2003. After his graduat ion he joined R&D department of AS- TR ON, the Netherlands Founda tion for Researc h in Astronomy , in Dwingeloo, T he Netherlands, w here he works with the system design and inte grati on group on the de ve lopment of the next gene ration of radio telesc opes. Since 2006 he is also with the Delft Uni ve rsity of T echnology , Delft , The Netherl ands, where he is pursuing a Ph.D. degre e. His research interests lie in the area of array signal processi ng, specifically calibr ation and imaging. VERSION 11 JUNE 2008 12 PLA CE PHO TO HERE Alle-J an van der V een (F’2005) was born in The Netherl ands in 1966. He recei ved the P h.D. degree (cum laude) from TU Delft in 1993. Throughout 1994, he was a postdoctoral scholar at Stanford Uni ve rsity . At present, he is a Full Professor in Signal Processing at TU Delft. He is the recipient of a 1994 and a 1997 IEEE Signal Processing Society (SPS) Y oung Author pa- per a ward, and was an Associa te Edito r for IEE E Tr . Signal Processing (1998–2001), chairman of IEEE SPS Signal P rocessing for Comm unicat ions T echni- cal Committee (2002-2004), and Editor-in-Chi ef of IEE E Signal Processing Letters (2002-2005). He currently is Editor- in-Chie f of IEEE Transa ctions on Signal Processing, and member-at-la rge of the Board of Gov ernors of IEEE SPS. His research interests are in the general area of system theory applied to signal processing, and in particula r algebra ic methods for array signal pro- cessing, with appli cati ons to wireless communicat ions and radio astronomy .