CORRELATION IMAGING Melvyn Wright, Lynn Urry, Matt Dexter, David - - PDF document

NSF, Washington, Aug 2008

CORRELATION IMAGING

Melvyn Wright, Lynn Urry, Matt Dexter, David MacMahon, Oren Milgrome, William Barott, Colby Gutierrez-Kraybill, Rick Forster, Garrett Keating, Gerald Harp, Dan Werthimer, Don Backer and the CASPER group Radio Astronomy Laboratory, Digital Group 1. SUMMARY

• what we have done.
• what we learned.
• what we plan to do.
• what we hope to learn.

– 2 – 2. CURRENT DIGITAL SYSTEM STATUS distributed data processing using multiple, re-usable hardware (1) ATA FX64 correlator – 100 MHz bandwidth, 1024 channels – custom F and X board design for up to 350 antennas – innovative, custom backplane solves signal distribution – pipeline correlator solves N 2 problem – 3 more FX64 correlators being built (Urry, MacMahon, Dexter) (2) CASPER – generic hardware and commercial protocols. – reuse hardware for IF processors, beamformers, correlators. – 10 Gb switch solves signal distribution (Werthimer) – flexible routing allows hardware to be maintained and upgraded.

• 42-antenna beamformer [100 MHz, + 1 × N correlator for calibration]
• 32-antenna beamformer [100 MHz, + 1 × N correlator for calibration]
• fly’s eye machine
• pulsar processor
• packetized correlator on EoR array (Backer & Bradley)

(3) SETI – analogue input from beamformer into PRELUDE. – 10 Gb digital input from beamformer into SONATA.

– 3 – 3. WHAT WE LEARNED

• (1) vindicate generic approach.

– custom designs take 5-10 years from concept to production.

• Multipurpose computing platform for radio telescopes

– system design philosophy treats boards as modular DSP resource. – faster to market. – allows user and student participation.

• FPGA enable wide range of radio astronomy applications.

– flexible interconnect architecture allows reconfiguration for multiple applications. – programming model allows focus on application rather than hardware.

• (2) beamformers calibration

– 1 × N correlator OK for L-band; not for X-band – need accurate calibration to form nulls.

• (3) RFI excision

– detailed editing can be automated. – big time saving for users. – enables real time imaging.

• (4) Data management

– sustained rate if we process at ATA and only keep “good” data. ∼ 100 GB/day = 10 hours processing on strato

– 4 –

• Fig. 1.— Bandwidth of Radio Astronomy correlators. Future correlators are shown in red.

The solid line shows Moore’s Law extrapolated from an SKA correlator for 4400 antennas with 1 GHz bandwidth. The right axis gives the data rate assuming 2000 spectral channels per GHz of bandwidth, 4 bytes per channel, a 25% overhead, and a 10s sample interval. ATA-42 ∼ 100 GB/day, ALMA ∼ 3 TB/day, SKA ∼ 1 PB/day

– 5 – 4. WHAT WE PLAN TO DO

• (1) Real Time Imaging
• Calibration and imaging integrated with hardware.

– stream data processing – flexible programming environment.

• Simultaneous imaging of multiple targets.
• Calibration in real time using global model.

– feedback calibration into beamformers and imagers.

• Subtract global model from uv data before imaging.
• Real time images update global model

– model becomes final image when observations completed.

• Calibrated images as normal output.
• (2) narrow-band software beamformer/correlator

– 64 ant × 2 pol × 20 MHz

• (3) wide-band gateware beamformer/correlator

– ATA correlator 64 ant × 2 pol × 600 MHz – CARMA correlator 23 ant × 2 pol × 8 GHz – VLBI beamformer 10 ant × 1 pol × 1 GHz

• (4) passive radar

– long delay and fast fringe rate in iBOB gateware.

– 6 –

• F

~

Packetized Control I F P

G

• F
• F

X X X

...

EthernetSwitch

I1

• )

ˆ ( ˆ ,

1

s V p f V

• Model(

) s,f,p PBeam( ) s,f,p Bandpass( ) f PolCal( ) f,p Gains( ) s,f,p

• Beam

Imager Astronomy Solver

I2

• )

ˆ ( ˆ ,

2

s V p f V

• k

j X Solver k j

• Fig. 2.— Data flow from telescopes to images.

– 7 – 5. WHAT WE HOPE TO LEARN

• Real Time Imaging

– transient astronomy – high fidelity imaging

• FPGA gateware survives by using a technology independent design flow.

– porting toolflow to new versions – porting gateware to new FPGA

• Best use of telescope and human resources.

– commensal observing

• How to build and operate SKA

– 8 –

• Fig. 3.— 10 pointings × 8 epochs of the ATATS field. These data were flagged using Garrett

Keating’s automatic editing. All 80 pointings were then mosaiced in MIRIAD to produce the image. The rms ∼ 4 mJy. Circles indicate NVSS sources down to 15 mJy, with the radius scaled with the NVSS flux density. Steve Croft 31July2008

– 9 – 6. MOTIVATION FOR REAL TIME IMAGING SCIENCE

• Delayed calibration and analysis limit science.
• Transient sources.
• Targets of opportunity.
• Observations directed to meet science goals.

TECHNICAL

• Non coplanar array geometry
• Non isoplanicity of atmosphere.
• Sources in sidelobes of beams.
• Beam pattern time variable
• RFI handled as data are acquired.
• Off-line data reduction limits data rates.

USER

• Telescopes used by non experts.
• Primary interest is astronomy, not data processing.
• Off-line data reduction expensive and time consuming.
• Expertize many astronomers do not have or want.
• Best use of telescope and human resources.

– 10 – 7. SYSTEM ARCHITECTURE Goal is to reduce data rate from sky to astonomer DIGITIZE

• Total data bandwidth is N × 2B × Npol × Nbits

4 1011 (N/100) (Npol/2) (B/GHz) (Nbits/8) bytes/s CHANNELIZE

• Large N, high dynamic range favors FX design.

– spectral resolution – bandwidth smearing – RFI mitigation (pre- and post correlation)

• Polyphase filter.

– Excellent channel separation. – Low cost ∼ log(Nchan) – Full bit growth then 4-bit selection for each channel.

• Parallel data processing.

– Rate in each channel by factor Nchan. e.g. with Nchan = 1000 2 108 (N/100) (Npol/2) (B/GHz) (Nbits/4) bytes/s/chan – calibration and imaging timing linear with channels and number of records. 10x as many channels == 10x as many records partitioning data over multiple processors is linear. up to disk I/O available; on strato is 600MB/s - 2GB/s

– 11 – CORRELATOR

• Data bandwidth for siderial source anywhere in sky:

107 (N/100)2 (Npol/4) (Nchan/103) (Nbytes/5) (Dmax/km) (λ/m)−1 bytes/s INTEGRATE AT MULTIPLE PHASE CENTERS

• Reduce data rate
• Range of fringe rates within the FoV limited by:

– primary beam, – isoplanatic patch, – non coplanar baselines, ∼ sqrt(λ/Dmax).

• Primary beam imaged by Nf ∼ λDmax/D2

ant images using 2D FFT.

e.g. ATA, with Dant = 6m, Dmax = 1 km, FoV defined by primary beam FWHM, ∼ 17 arcmin at λ = 3 cm, and by non coplanar baselines at λ = 1 m.

• CMAC for each baseline and frequency channel, but data rate reduced.

– data bandwidth for imaging primary beam FWHM: 106 (N/100)2 (Npol/4) (Nchan/103) (Nbytes/5) (Dmax/km) (Dant/m)−1 bytes/s

• REAL TIME IMAGING

– 1 image per hour

– 12 –

• Fig. 4.— Multichannel Calibration and Imaging

– 13 – 8. CALIBRATION

• Calibrate using global sky model.
• Separately calibrate each isoplanatic patch.
• Global calibration model across array.
• Image regions of interest and confusing sources.
• Stream data processing.

Algorithm

• Calculate model visibility from sky model

V

′

j,k = exp(2πi/λ r.s0) × Σ I × A × B × P × G × exp[2πi/λ r.(s − so)]

– I(s, ν, p) sky model. – A(s, ν, p) primary beam, – B(ν) instrumental bandpass, – P(s, ν, p) polarization calibration, – G(time, s0) gain versus time and phase center. – r antenna baseline vector. – s, ν, and p position, frequency and polarization. – so is the phase center for each region of interest.

• Measurement equation for A, B, P, and G.

– decompose into antenna dependent components Gjk = gj × gk. – use a-priori measurements wherever possible.

• Global sky model to determine antenna gains g(t, s0).

– gains measure tropospheric and ionospheric delays.

• Self calibration. χ2 = < Σ[V × gigj − V ′]2 / σ2

ij >

– χ2 accumulated in distributed processors.

– 14 – 9. IMAGING

• Parallel processing in distributed architecture.
• Subtract a-priori source model from calibrated uv data.

– confusion and sidelobe subtraction.

• Image each phase center.

– image size < Dmax/λ for 2D FFT, ∼ 108 for 1000 km baseline at λ 1 m.

• Deconvolution minimized by good uv coverage for large N.
• Imaging guided in real time by convergence of model and χ2 image.

– variable sources are inconsistent with the global model. – χ2 image identifies transient sources.

• Phase centers can be moved to image regions for science goals,

– image confusing sources. – image different isoplanatic patches.

• Images update global model.
• Model and calibration improve as observations proceed.

– observe until model consistent with uv data streams.

– 15 – 10. APERTURE ARRAYS

• Antenna size, Resolution and Confusion.
• A single dish collecting area ∼ ηD2; resolution FWHM ∼ λ/D.
• Large antennas are confusion limited.

– enough sensitivity to detect more than one source within the beam. = enough collecting area, but need more resolution => build an array: collecting area ∼ ηND2; resolution FWHM ∼ λ/dmax.

• Need to add up the E-field across distributed aperture - preserving the relative phase
• f the wavefront across an array of antennas.

– This is quite difficult. – electronics at each antenna need to preserve phase, – atmospheric phase shifts distort the wavefront, – the problem is rather like adaptive optics. – keep the path lengths within ∼ λ/20 to make an accurate telescope.

– 16 – 11. BEAMFORMERS and CORRELATORS

• Array aperture is set of antennas each with its own E-field and phase.

Beamformer Voltage response = Σ Vj exp 2πi λ bj.s ds Beamformer Power = Σ V 2

j + Σ Vj V ∗ k

exp 2πi λ (bj − bk).s ds sum of total powers from all antennas + cross products of antenna pairs Beamformers

• Phased array beams into expensive backends.
• Multiple compact targets ( pulsars, SETI,... )

=> need a correlator to calibrate beamformer. Correlators

• Sampling < Vj V ∗

k > allows imaging full FOV of antennas.

• Image multiple FOVs simultaneously.
• 40 year development of calibration and imaging techniques.

– signal distribution => custom backplane or switch. – high data rate => real time data processing.

• 4000 antenna, 1 GHz correlator feasible for SKA by 2020.

– 17 – 12. QUESTIONS, EXTRA SLIDES ETC. 13. Acknowledgments Thanks to the dedicated efforts of many people who make this possible, in particular: Lynn Urry, Matt Dexter, David MacMahon, Oren Milgrome, William Barott, Colby Gutierrez-Kraybill, Rick Forster, Garrett Keating, Gerald Harp. Dan Werthimer, Don Backer and the CASPER group. Students at summer schools and imaging seminars for their enthusiasm and energy.

– 18 – 14. APERTURE SYNTHESIS IMAGING

• Array response

V (t) =

• I(s)A(s − s′) exp 2πi

λ b.s ds I(s) is the sky brightness distribution A(s − s′) is the primary beam illumination of the sky s′ is pointing center

• Phase tracking array

V = exp 2πi λ b.s0(t)

• I(σ)A(s − s′) exp 2πi

λ b.σdσ Instrumental terms Source structure s0 is phase tracking center σ = s − s0 is the vector from the phase tracking center to the source σ = (x, y, z); b = (u, v, w) – an interferometer array is a chromatic instrument – use a small range of λ or else we get bandwidth smearing. => continuum sources can use bandwidth synthesis (MFS) to obtain more Fourier samples of the brightness distribution.

– 19 –

• Imaging
• We only have discrete samples of V, so we define a weighting function W, and evaluate:

I′(x, y) =

• W(u, v)V (u, v) exp −2πi

λ (ux + vy) dxdy – weighting function W can have any value where V is sampled, and W = 0 if not. – I′ is Fourier transform of product of W and V , = convolution of the Fourier transforms of W and V : I′(x, y) = B(x, y) ⋆ [I(x, y)A(x, y)] where B(x, y) =

• W(u, v) exp −2πi

λ (ux + vy) dudv B is the synthesized beam.

– 20 – 15. SIDELOBE SUBTRACTION Problem

• Sidelobes of sources outside the regions of interest confuse the images,

– full FoV cost ∼ λ1.7NstaD3

max/D6 ant, Nsta = antennas/station. (Perley & Clark).

– non-coplanar baselines, ∼ λN 2D3

max/D4 ant, (Cornwell)

– for constant collecting area, ∼ λD3

max/D8 ant.

• These scaling relations only apply to full FoV at highest resolution.

(Lonsdale, Doeleman, & Oberoi, 2004; Wright, 2004).

• Larger antennas do not solve the problem.

– poor uv coverage and mosaicing degrade image fidelity.

• Off-line processing expensive for large N arrays.
• RFI

Solution

• Subtract global model from uv data stream.

– stream processing of highly parallel data. – distributed processors for each region of interest.

• Calibration includes primary beam, bandpass, polarization.
• Characterize RFI as function of time, frequency and polarization.
• RFI subtracted from uv data stream before beamformers.

– 21 – 16. ARCHIVE

• uv data streams and metadata saved in archive.
• uv data streams from archive replayed though imaging system
• improve calibration using best model of sky and phase screen.

UVDATA

• ATA64 x 1 pol x 4 bytes/sample x 1024 chan x 1s sample x 3600 x 24 hours x 4 FX64’s
• Survey Science - wide field/no fringe tracking. 1 s dump

– 2.8 TB/day – dump correlator at 1s rate, do flagging and RFI excision on 1s data – integrate to 10s in multiple FOV (== post correlation fringe tracking – 280 GB/day/phase center

• non-survey science targets

– primary beam FOV, Nyquist sampling. 28s dump ∼ 100 GB/day – modest FOV with fringe tracking. 10s dump ∼ 280 GB/day

– 22 – IMAGES

• Real time imaging

– keep the original uv-data for limited time. – until final, best calibration has been done. – until the image data products have been verified.

• Survey Science

– image data rate = 20482 x 4 byte x 24 hours x 4 pols x 4 bands – ∼ 6.4 GB/day

• non-survey science targets

– depends on the science; ∼ 1 image/hour. – rate is about the same as survey science

• Spectral Line Images.

– 10-20 % of observations need to keep the spectral info, – ∼ 200 good channels, then images are 200x bigger. – ∼ 128 GB day for spectral line images.

• BOTTOM LINE

∼ 100 GB/day sustained rate if we process at ATA and only keep “good” data.

– 23 – PROCESSING SPEEDS SunOS strato i86pc i386 i86pc Solaris

• UVCAL (read+write Miriad format )

66 MB/s

• calibration:

UVGEN compute model + Miriad write to disk 2.8 MB/s

• multichannel imaging

INVERT: Miriad read, write scratch file, read scratch file, FFT, write image 7.0 MB/s

• disk copy

cp -r ata-42.30.uv junk about 1s. on strato (/ataarchive/scratch/obs/mchw/ata) = 829 MB/s 9s on sun (sun1 disk) = 92 MB/s 70s on boss ( /hp disk) = 12 MB/s

• BOTTOM LINE

calibration + 2.8, imaging 7.0 MB/s ∼ 10 GB/hour 100 GB/day = 10 hours processing on strato

– 24 – 17. PROOF OF CONCEPT

• Radical departure from conventional off-line data reduction.

(1) Simulated data

• Image uv data at multiple phase centers.
• Calibrate at each phase center.
• Subtract model at each phase center.
• Image each phase center.
• Simulations show we can deconvolve multiple regions from uv data.

– source subtraction degrades slowly in FoV ∼ sqrt(λ/Dmax). – offline simulation multiple regions are handled sequentially.

• Goal is to develop efficient handling of uv data streams.

– use existing software in floating point processors.

• Distributed processors for each phase center.

– 10-20,000 uv samples/s (1-3 Mbytes/s) on modest sized desktop

– 25 –

• Fig. 5.— A mosaic MFS Image of CASA scaled up 40 times real size, imaged with the ATA

at 1.42 GHz. The image is 4 degrees on a side. The synthesized beam FWHM, 77 x 78 arcsec is shown in the lower left corner.

– 26 – (2) Real Data

• EPOCH OF REIONIZATION EXPERIMENT

(Backer and Bradley, 2005)

• Dipole array 150-200 MHz, 513 channels sampled at 0.1 Hz,
• uv data phase rotated to Cygnus A and Cas A

– antenna positions fitted over a limited HA range.

• Self-calibrate at each phase center.

– no correction for primary beam response or polarization

• Dynamic range ∼ 100
• Spectral dynamic range ∼ 1000

– 27 –

• Fig. 6.— Cygnus Image at 175 MHz from 4-dipole array

– 28 –

• Fig. 7.— Cas A Image at 175 MHz from 4-dipole array

– 29 –

• Fig. 8.— FX64 correlator block diagram

– 30 –

• !

" " #$ " %" &'" " (

• $$

)"! '"

• * !

" " * *

• $"

&+"",&- #,)." /,"(."""&,"&'." '""))&+""#$) $$ &+""0,) /,")"!("&&-"0,) /,"'"0,) 01("+"!'".",-*) * !&""*) ** !&""*))02) #,""1.3-"1&'"&"""*) 0*&'" &" " " &"& * *,!&'"&"""&"&*+"4)

• Fig. 9.— FX64 correlator block diagram

– 31 – 18. STREAM DATA PROCESSING

• data processing is integrated with DSP hardware
• Block floating used in hardware DSP.
• MIRIAD uv data is a stream of scaled 16-bit integers.

– metadata is a stream of named variables and values. 19. HARDWARE PROTOTYPE

• Multipurpose computing platform for radio telescopes
• FPGA enable wide range of radio astronomy applications.

– BEE2 250 109 CMAC/s, memory 12 GB/s, I/O 360 Gb/s

• Custom designs take 5-10 years from concept to production.
• System design philosophy treats boards as modular DSP resource.

– flexible interconnect architecture allows reconfiguration for multiple applications. – programming model allows focus on application rather than the hardware. – software survives by using a technology independent design flow.

• Programmable digital lens for radio astronomy.

– 32 –

• Fig. 10.— Beamformer block diagram.

– 33 – 20. SUMMARY

• Calibration and imaging integrated with hardware.
• Stream data processing
• Flexible programming environment.
• Simultaneous imaging of multiple targets.
• Calibration in real time using global model.
• Calibration feedback into beamformers and imagers.
• Subtract global model from uv data before imaging.
• Real time images update global model
• Model becomes final image when observations completed.
• Calibrated images as normal output.
• Best use of telescope and human resources.

– 34 – Table 1: FOV AND SAMPLE RATES FOR ATA AND SKA sqrt(λ/Dmax) is the limit to the image size imposed by non-coplanar baselines. Nf ∼ λDmax/D2

ant is the number of 2D FFT

c/Dmax is max channel width without delay tracking. Fringe is the fringe rate. Nyquist is the sample rate within primary FWHM. λ Dmax Dant λ/Dmax λ/Dant sqrt(λ/Dmax) Nf c/Dmax Fringe Nyquist [m] [km] [m] [arcsec] [arcmin] [arcmin] [kHz] [Hz] [Hz] 1.00 1 6 206.26 572.9 108.7 27.8 300 0.07 0.02 0.21 1 6 43.32 120.3 49.8 5.8 300 0.35 0.02 0.03 1 6 6.19 17.2 18.8 0.8 300 2.42 0.02 1.00 1 12 206.26 286.5 108.7 6.9 300 0.07 0.01 0.10 1 12 20.63 28.6 34.4 0.7 300 0.73 0.01 0.01 1 12 2.06 2.9 10.9 0.1 300 7.27 0.01 1.00 10 12 20.63 286.5 34.4 69.4 30 0.73 0.12 0.10 10 12 2.06 28.6 10.9 6.9 30 7.27 0.12 0.01 10 12 0.21 2.9 3.4 0.7 30 72.7 0.12 1.00 100 12 2.06 286.5 10.9 694.4 3 7.27 1.21 0.10 100 12 0.21 28.6 3.4 69.4 3 72.7 1.21 0.01 100 12 0.02 2.9 1.1 6.9 3 727. 1.21 1.00 1000 12 0.206 286.5 3.4 6944. 0.3 72.7 12.12 0.10 1000 12 0.021 28.6 1.1 694.4 0.3 727. 12.12 0.01 1000 12 0.002 2.9 0.3 69.4 0.3 7270. 12.12 1.00 1000 25 0.206 137.5 3.4 1600. 0.3 72.7 5.82 0.10 1000 25 0.021 13.8 1.1 160. 0.3 727. 5.82 0.01 1000 25 0.002 1.4 0.3 16. 0.3 7270. 5.82

– 35 – 21. References Backer and Bradley, 2005, NRAO newsletter January 2005, no. 102, p27. Brodersen, B., Chang, C., Wawrznek, J., Werthimer, D., & Wright, M., 2004 ”BEE2: A Multi-Purpose Computing Platform for Radio Telescope Signal Processing Applications” http://bwrc.eecs.berkeley.edu/Research/BEE/BEE2/presentations/BEE2 ska2004 poster.pdf Cornwell, T.J. & Perley, R.A., 1992, ”Radio-Interferometric Imaging of Very Large Fields”, A&A 261, 353 Cornwell, T.J., Holdaway, M.A. & Uson, J.M., 1993, A&A 271, 697, ”Radio-interferometric imaging of very large objects: implications for array design”, Cornwell, T.J., Golap, K., & Bhatnagar, S., 2003, ”W projection: a new algorithm for non-coplanar baselines”, EVLA memo 67 Cornwell, T.J., 2004, ”EVLA and SKA computing costs for wide field imaging (Re- vised)” EVLA memo 77 Jones, D.L., 2003, ”SKA Science Requirements”, SKA memo 45, version 6, (D.L.Jones, 16 December 2003, ISSC11-9.4) Lonsdale, C.J., Doeleman, S.S., & Oberoi, D., 2004, ”Imaging Strategies and Postpro- cessing Computing Costs for Large-N SKA Designs”, SKA memo 54, July 2004. Perley, R. & Clark, B., 2003, ”Scaling Relations for Interferometric Post-Processing”, EVLA memo 63 Thompson, A. R., Moran, J. M. & Swenson, G. W., 2001, ”Interferometry and synthesis in radio astronomy”, 2nd ed. New York : Wiley, 2001. Wright, M.C.H., 2002, ”A model for the SKA”, SKA memo 16 Wright, M.C.H., 2002b, ”Allen Telescope Array Imaging”, ATA memo 52. Wright, M.C.H., 2004, ”SKA Imaging”, SKA memo 46 Real Time Imaging, Melvyn Wright, 2005, SKA MEMO 60 SKA Task Force on Survey optimization, Melvyn Wright(chair), 2006, SKA MEMO 81 Surveys, Calibration and Imaging, Melvyn Wright Invited talk, 2006, SKA Joint Work- ing Group. Image Fidelity as a Function of Source Size and Calibration Errors, Melvyn Wright,

– 36 – 2007, CARMA memo 38 Deconvolving Primary Beam Patterns from Mosaic and Polarization Images, Melvyn Wright & Stuartt Corder, 2008, CARMA memo 43