Import standard modules:


In [ ]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import HTML 
HTML('../style/course.css') #apply general CSS

Import section specific modules:


In [ ]:
from IPython.display import Image
import scipy.signal

In [ ]:
HTML('../style/code_toggle.html')

7.3 Analogue Electronics (G- and B-Jones)

In the first sections of this chapter we have introduced Jones notation and the radio interferometric measurement equation (RIME). These ideas have been presented in the abstract as a way to lay out the mathematics of how we describe observing systems. For the remainder of the chapter we will explore the details of different instrumentational effects (and the corresponding Jones matricies), and how these effects alter the sky signal. The recovering of the original sky signal, which is the ultimate goal, requires knowledge of the observation system. How well we understand an observing system affects how accurate we can recover the true sky signal. We will discuss the process of calibration in the next chapter.

In $\S$ 7.2.3 ➞ we introduced the concepts of direction-dependent and direction-independent effects (DDE's and DIE's). Historically, direction-independent effects were the primary concern during observations. This was mainly due to simplicity, hardware and array limitations, and because DIE's are the dominate effects on the sky signal. These effects occur at the point on conversion of the sky electric field to a voltage at the feed and propogates through the telescope and array electronics. Typically, all the DIE's are wrapped up into a single 'gain' Jones matrix called the $\mathbf{G}$-Jones. This term is made up of numerous effect, primarily due to the analogue electronic components such as filters, amplifiers, mixers, cables, etc. This will be discussed through the section.

Often, there is a second DIE Jones matrix, the bandpass or $\mathbf{B}$-Jones which is related to the $\mathbf{G}$-Jones. As we have discussed in the $\S$ 7.1 ➞ a Jones matrix is a function of time and frequency. But, it is often a good approximation to think of the analogue electronics effects to be separable into a frequency-variable, time-stable function (the bandpass) and a time-variable, frequency-stable function (gain), such that the total, time- and frequency-dependent gain matrix $\mathbf{G}'$ is

$$ \mathbf{G}'(t, \nu) \approx \mathbf{G}(t) \cdot \mathbf{B}(\nu) $$

Thus, confusingly, the literature and software the $\mathbf{G}$ can be either a time-dependent or time- and frequency-dependent term depending on how the author defined it. Be sure to be aware of this inconsistency. The most important point to make clear at this moment is that the $\mathbf{G}$ term has the largest effect on the observed signal compared to other Jones terms, which is why we are starting with it first.

7.3.1 An Idealized Source

For the moment we will ignore interferometry and consider only a single radio telescope. Consider a simple sky which contains a single source, what would we hope to observe? The source would have some intrinsic flux $I$ (measured in Janskys) and spectral index $\alpha$ of the form

$$ I(\nu) = I_0 \left( \frac{\nu}{\nu_0} \right)^{-\alpha} $$

where $I_0$ is the flux of the source at some reference frequency $\nu_0$. We think of most radio sources as being stable over long lifetimes, thus the source flux will not noticably change over short periods of time. Though the signal we are measuring is a time-domain signal, we will look at the frequency spectrum through out this section as that is what we are ultimately interested in. If we were to plot the flux of the source as a function of time and observing frequency in a waterfall plot it would look like the figure below


In [ ]:
# Setup an ideal source coherency spectrum
nchans = 512 # number of frequency channels
freq0 = 1.1e9 # start frequency (Hz) (centre of first channel)
freq1 = 1.5e9 # stop frequency (Hz) (centre of last channel)
freqs = np.linspace(freq0, freq1, nchans)
ntime = 500

I0 = 1. # flux at reference frequency (Jy)
nu0 = 1.421e9 # reference frequency (Jy)
alpha = -0.25 # source spectral index
Inu = lambda nu: I0 * (nu/nu0)**alpha # Stokes I spectrum

# Stokes Parameters
I = Inu(freqs)
Q = np.zeros_like(I)
U = np.zeros_like(I)
V = np.zeros_like(I)

# Compute the polarisation values
brightness = np.empty(shape=(nchans, 2, 2), dtype=np.complex128)
brightness[:,0,0] = I + Q
brightness[:,0,1] = U + V*1j
brightness[:,1,0] = U - V*1j
brightness[:,1,1] = I - Q

brightness = np.tile(brightness[np.newaxis], [ntime, 1, 1, 1])

In [ ]:
fig, axes = plt.subplots(figsize=(8,8))

ax1 = plt.subplot2grid((5,5), (0,0), colspan=4, rowspan=4)
plt.imshow(np.abs(brightness[:,:,0,0]+brightness[:,:,1,1]), aspect='auto')
ax1.get_xaxis().set_visible(False)
ax1.get_yaxis().set_visible(False)

ax2 = plt.subplot2grid((5,5), (4,0), colspan=4, rowspan=1)
plt.plot(freqs/1e9, brightness[10,:,0,0]+brightness[10,:,1,1])
plt.ylabel('flux (Jy)')
plt.xlabel('Frequency (GHz)')

ax3 = plt.subplot2grid((5,5), (0,4), colspan=1, rowspan=4)
plt.plot(brightness[:,10,0,0]+brightness[:,10,1,1], np.arange(brightness.shape[0]))
ax3.get_xaxis().set_visible(False)
ax3.invert_yaxis()
ax3.yaxis.tick_right()
plt.ylabel('Time (s)')
ax3.yaxis.set_label_position('right')

plt.suptitle('Idealized Source Spectrum')

The bottom plot shows the spectrum for a single moment in time, and the plot on the right shows the flux at a single frequency. The source is stable in time and frequency. This is what we would like to detect. But, of course, things are not so simple. Any system we build will introduce noise into the measurement. And, as we have hinted at, our measurement systems alters the original signal. Instead of this nice, ideal measurement we see something like the waterfall plot at the end of this section. We will step through the various stages of an example analogue front-end to understand why the observed spectrum differs from the original source spectrum.

7.3.2 System Noise and Sensitivity

All measurement systems have noise. In fact, anything which is not at absolute zero (that is everything) has some energy, which is radiated as heat, another term for heat is noise. We are interested in measuring some signal. But, if our system introduces too much noise into the measurement then we can not recover the original signal. When we discuss the sensitivity of a telescope we are primarily interested in how much noise is added from the system during the measurement. This system noise is a result of, anong other issues, electronic components, the engineering design of the telescope (antenna feed, dish design, etc.), and sky signals. Newer electronic components tend to have lower noise characteristics compared to older systems. Thus, telescopes are often upgraded with newer systems.

Noise as Temperature

It is very common to discuss noise as a temperature. As we stated, noise is another term for heat, and heat is measured as a temperature. From an electronics point of view we can consider a signal to be a voltage in an electrical circuit. The signal is then a Johnson-Nyquist noise source, the power dissapated by that noise source is

$$P = k_B T \Delta \nu$$

Where $P$ is the power of the voltage signal, $k_B$ is the Boltzmann constant, $T$ is temperature in Kelvin, and $\Delta \nu$ is the bandwidth of the signal in Hertz. Thus, in literature the noise power of a radio telescope (or a component of one) is often described in terms of a temperature

$$T = \frac{P}{k_B \Delta \nu}$$

System Temperature

Given that we would like to determine the sensitivity of a telescope we first need to know the noise characteristics. Then, we can decide if it is possible to detect the signal we are after. We start by saying there is some temperature which defines the noise characteristics of the telescope, as one would expect this is the system temperature and is made up of multiple components

$$ T_{\textrm{sys}} = T_{\textrm{sky}} + T_{\textrm{atmosphere}} + T_{\textrm{spillover}} + T_{\textrm{rx}} + \ldots$$

where $T_{\textrm{sky}}$ is the stocastic signal from the sky. At L-band (1.4 GHz) the sky temperature is approximately 10 K, 2.73 K of which are from the Cosimic Microwave Background (CMB). The temperature of the sky varies depending on the observing frequency, for example at low frequency, say 100 MHz, the sky is close to 1000 K primarily due to synchrotron radiation. We can approximate the 'room temperature' of Earth as 300 K. So, if the sky noise is significantly 'hotter' than room temperature then there is little advantage to building instruments with low noise characteristics.

There is also noise, again depending on observing frequency, due to the atmosphere and ionosphere $T_{\textrm{atmosphere}}$. Though this noise is naturally produced we still think of it as being part of the observing system. If our telescope was placed on the moon then there would be essentially no atmospheric noise. We discuss further details of different propagation effects in $\S$ 7.7 ➞. At L-band this is a small noise effect ($\sim 1$ K), the phase rotation effects discussed later in the chapter are a much more serious effect.

The Earth is a 300 K black-body radiation source, i.e. a very hot noise source. We would like to block the noise from the earth (called ground noise) entering our measurement system. This can be done a number of ways, such as cooling the electronics, and isolating the system. One way the ground noise enters the system is from the illumination pattern of the primary beam. We discuss primary beams later in the $\S$ 7.5 ➞. But, for now it is important to know that the telescope optics produce a 'beam' which is sensitive to all locations on the sky. The optics are designed such that the ebam is most sensitive to the direction the dish is pointed, but there is always some sensitivity in other directions. Thus, no matter where in the sky the telescope is pointed it is sensitive to the picking up signal from the ground. If the telescope is pointed low on the horizon, an even more sensitive section of the ebam will be pointed at the ground. This is called spillover noise $T_{\textrm{spillover}}$ and is a pointing direction-dependent noise term. If the integrated spillover is on the order of a few precent then this can result in noise of 10's of Kelvin.

The final system temperature term of note is the receiver noise $T_{\textrm{rx}}$, this is the noise introduced inot the measurement by the analogue electronics used to convert the weak sky signal into a digital signal which cna be recorded. A KAT-7 or MeerKAT telescope produce approximately 20 K in receiver noise. We will discuss the analogue receiver system further in the rest of this section.

Radiometer Equation

For a single radio telescope the measurement of the average power over some bandwidth for a direction in the sky is done with a radiometer. Given that the noise signal we are measuring is approximately Gaussian, the error of the measured signal can be described by the ideal radiometer equation

$$ \sigma_{T} = \frac{T_{\textrm{sys}}}{\sqrt{\Delta\nu \tau}} $$

where $\sigma_{T}$ is the standard deviation of the measured noise temperature about the true noise temperature. The radiometer equation describes the necessary bandwidth ($\Delta \nu$) and integration time ($\tau$) required to reach a desired noise level ($\sigma_{T}$) for an ideal system with a given system temperature ($T_{\textrm{sys}}$) when observing a broadband source.

From the radiometer equation we can see that a linear decrease in system temperature will lead to a linear reduction in $\sigma_{T}$. A linear increase in the bandwidth or intergration time only leads to square-root improvement in $\sigma_{T}$. For example, If there are two systems, with the same bandwidth $\Delta \nu$, in which one has a system temperature twice that of the other, i.e. $T_{\textrm{sys,0}} = 2 T_{\textrm{sys,1}}$, then $T_{\textrm{sys,0}}$ will take 4 times longer to reach the same sensitivity as $T_{\textrm{sys,1}}$.

The radiometer equation has been presented for an ideal system, but in reality analogue systems are unstable and introduce 'gain fluctuations'. These fluctuations come about due to a number of effects: the physical temperature of the electronics, cross-coupling leakage, RF sheilding among others. The stability of the system during an observation will depend on how well these effects are managed. These gain fluctuations will have both a time and frequency dependence.

System Equivalent Flux Density (SEFD)

The overall sensitivity of a radio telescope is described by the system equivalent flux density (SEFD) which is a ratio of the noisiness of the instrument ($T_{\textrm{sys}}$) to the effective signal gain $G_{eff}$. Where gain is a term for all the various amplifications and attenuations of the sky signal we are interested in. Formally the SEFD is defined as

$$ \textrm{SEFD} = \frac{T_{\textrm{sys}}}{G_{eff}} = \frac{2 k_B \eta T_{\textrm{sys}}}{A_{eff}} $$

where $k_B$ is the Boltzman constant, $\eta$ is an generic efficiency factor (always greater than 1) used to represent the signal loss in various components of the telescope (e.g. dish surface inefficiencies, feed under-illumination, quantization, etc...), and $A_{eff}$ is the effective collecting area of the telescope. For a dish telescope $A_{eff}$ can, to first order, be approximated as the geometric surface area of the dish.

A decrease in the SEFD is equivalent to an improvement in the sensitivity of the telescope.

7.3.3 Analogue Receiver Front-End

To understand the details of the system temperture and gain we need to start with the analogue receiver front-end, also called the analogue chain. This is the set of analogue electronic components which are used to condition the weak sky signal in order to do digital computations (like correlation). Ideally we would like to have a signal in a digital form for a number of reasons. Digital signals are predictable, and any computation is deterministic. It is easier to do computation like Fourier transforms and correlations on digital signals compared to analogue signals. Analogue signals have the advantage that they are continuous, there are no Nyquist sampling or aliasing issues. But, analogue components are effected by feedback, cross-talk, temperature, and power supplies. And, worst of all, analogue components introduce noise into the signal. In general, analogue components produce variable response in time and frequency where as a digital system has a stable response in time and frequency.

We would prefer to digitize our sky signal as soon as possible. Unfortunately, the sky signal is very weak and needs to be amplified and conditioned before it can be converted to digital. The analogue-to-digital converter (ADC) transforms an analogue signal to a digital signal, but if the signal is too weak then the ADC can not sample the analogue signal. Below is a block diagram of an example analogue front-end.

Amplification is the effect of increasing the strength of a signal absolutely or relative to another signal. In electronics an amplifier might be used for this, but there are other types of amplification such as building a larger dish to collect more signal. The opposite effect is called attenuation. All Jones matrices can be thought of performing an amplification or attenuation on the input signal to produce the output signal. This amplification (attenuation can also be seen as an amplification less than 1) can be complex or real-valued, a function of time, frequency, or direction. Some Jones matrices will have unity gain such as a rotation transform. Signal conditioning is the application of different filters to select out a limited signal of interest.

Figure: Block diagram example of an analogue receiver front-end for a radio telescope. A current is induced in the receptor by the electric field of the sky source. The resistance of the components creates a voltage which is propotional to the amplitude of the electric field. This voltage is amplified by a low noise amplifier (LNA) in a cryostat at 20 K. A second stage amplifier further amplifies the signal in a 70 K intermediate cryostat. The signal is bandpass filtered to select the band of interest. A bandstop filter is used to filter out frequencies with radio frequency interference (RFI). A local oscillator (LO) is mixed with the signal to down mix the original radio frequency (RF) to an intermediate frequency (IF) which is then digitized by an analogue-to-digital converter (ADC).

The electric field of the sky signal induces a current on the receptor (also known as a feed),we will discuss feeds in $\S$ 7.4 ➞. This is a very weak signal, we need to use an amplifier to increase the strength of the signal. The process of amplification introduces noise, we would like to minimize this noise. As technology developments improve the noise introduced by amplifiers has decreased. The common method to reduce the noise of an amplifier is by cooling it using a cryostat. This is a chamber under vacuum which can be cooled down to very low temepratures (20 K if using liquid helium, 70 K for liquid nitrogen). Remember, noise and temperature are the same, so the lower the temperature the lower the noise. This special amplifier is called a low noise amplifier (LNA) and is usually the main component of the system noise. When we are building a telescope system we want to focus on making sure the LNA has the best noise perforance possible (i.e. the LNA is a significant cost to the analogue front end).

The idea of the LNA is that we have a weak signal which needs to be amplified with the smallest amount of noise added into the signal. Once the signal has been amplified with the LNA we can use normal amplifier and apply filters because the signal to noise ratio will be sufficiently large that the added noise from these components will not greatly affect the original signal. This can be seen in the receiver temperature equation.

$$ T_{\textrm{rx}} = T_{\textrm{feed}} + \frac{T_{\textrm{passive}}}{G_{\textrm{feed}}} + \frac{T_{\textrm{LNA}}}{G_{\textrm{feed}} G_{\textrm{passive}}} + \frac{T_{\textrm{amp}}}{G_{\textrm{feed}} G_{\textrm{passive}} G_{\textrm{LNA}}} + \ldots$$

where the $T$ terms are the additional noise temperatures introduced with each component, and $G$ is the 'gain' term which can be less than 1 (attenuation) or greater than 1 (amplification). The gain from the feed and passive components will be less than 1, and are similar to an efficiency term in when discussing the primary beam, these will cause the temperature from these components to increase in the system temperature. By having the LNA as early as possible we see that the weak sky signal is amplified, and noisier components such as filters and second stage amplifiers can be introduced without introducing significant noise compared to the sky signal.

From the receiver temperature equation above we can see that each additional temperature term includes the gain of the LNA in the denominator. Thus, the more gain introduced by the LNA the smaller each additional noise term. There are some noise terms which come before the LNA stage. There is always a temperture of the feed $T_{\textrm{feed}}$ (any physical feed will have some resistance or else there would be no current). There are passive components such as cables, orthomode transducers (OMTs), waveguides which all add a noise component $T_{\textrm{passive}}$.

Decibel Units

As a quick aside, it is common to discuss electronic and system gain using logarithmic decibel units. This is worth a quick review. We usually describe a signal in linear units. But, it is often more intuitive to use logarithmic units when we are discussing values which cover many orders of magnitude such as in electronics. The standard unit to use is decibels which for a power $P$ is defined as

$$ P_{dB} = 10 \log_{10} \left ( \frac{P}{P_0} \right ) $$

where $P_0$ is the reference power, typically set to unity. Important values to remember are: 3 dB is a factor of 2 increase, 10 dB is a factor of 10, -3 dB is a factor of 0.5. Every factor of 10 increase in decibels is an order of magnitude increase. Note, if the signal is not a power, such as a voltage or electric field the definition of decibel is

$$ P_{dB} = 20 \log_{10} \left ( \frac{V}{V_0} \right ) $$

as we need to compute the power of the voltage in the process of converting to decibels. This is useful for describing the response of Jones matricies in decibels.

Back to the LNA in our analogue front-end example, it is typical to amplify the sky power signal by 20 to 30 decibels with the LNA. That is, the LNA amplifies the sky signal by a factor of 100 to 1000. Now, we can use other components and not need to worry about introducing significant noise relative to the signal strength.

Below is an example of the LNA response in time and frequency. This is a simple model where the LNA bandpass response is stable in time, but the overall system gain varies in time. The bandpass has a response that ranges between 30 dB and 20 dB depending on the frequency. There is a'ripple' across the band due to the LNA design. Overall the LNA has a peak around 1.2 GHz and drops off at higher frequency. The overall gain varies a few precent over time, this can be due to a number of effects. The heating and cooling of the cryostat is never perfect so there is variation the gain. Depending on the observing source we should expect the gain to change, a bright source will introduce more noise than a dimmer source. Feedback and cross-talk between analogue components will always exist even with good isolation. This is a farily typical LNA response one would expect to see.


In [ ]:
#LNA bandpass
#xrange: 0,1
#yrange: 20-30 dB
#a made up LNA response spectrum (in dB)
nfreqVals = np.array([-0.5, 0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
                  0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.05])
lnaVals = np.array([27, 27, 28.1, 28, 27.5, 29, 30.1, 30, 28.5, 28.5, 30,
                  27, 27, 25.5, 23, 24, 25, 22.5, 21, 21, 22, 19.5, 19.5])

#add a little noise for randomness
lnaStd = 0.5 #standard deviation of noise in dB
noiseVals = np.random.normal(loc=0.0, scale=lnaStd, size=lnaVals.shape)

nfreqs = np.linspace(0, 1, num=nchans) #normalized frequencies
coeffs = np.polyfit(nfreqVals, lnaVals + noiseVals, deg=15)
lnaGain = 10.**(np.polyval(coeffs, nfreqs)/10) #convert to linear

#Gain variation
#xrange: 0,1
#yrange: 1 + std
tVals = np.linspace(-0.5, 1.5, num=60)
gainStd = 0.10
gVals = np.random.normal(loc=1.0, scale=gainStd, size=tVals.shape)

ts = np.linspace(0, 1, num=ntime)
gCoeffs = np.polyfit(tVals, gVals, deg=20)
gain = np.polyval(gCoeffs, ts)

gainSpectrum = np.outer(gain, 10. * np.log10(lnaGain))

fig, axes = plt.subplots(figsize=(8,8))

ax1 = plt.subplot2grid((5,5), (0,0), colspan=4, rowspan=4)
plt.imshow(np.abs(gainSpectrum), aspect='auto')
ax1.get_xaxis().set_visible(False)
ax1.get_yaxis().set_visible(False)

ax2 = plt.subplot2grid((5,5), (4,0), colspan=4, rowspan=1)
plt.plot(freqs/1e9, gainSpectrum[100, :])
plt.ylabel('Gain (dB)')
plt.xlabel('Frequency (GHz)')

ax3 = plt.subplot2grid((5,5), (0,4), colspan=1, rowspan=4)
plt.plot(gainSpectrum[:, 100], np.arange(gainSpectrum.shape[0]))
ax3.get_xaxis().set_visible(False)
ax3.invert_yaxis()
ax3.yaxis.tick_right()
plt.ylabel('Time (s)')
ax3.yaxis.set_label_position('right')

plt.suptitle('LNA Response')

Filters

It is necessary to apply signal filters to the observed signal for a number of reasons. We are limited to a finite observing bandwdith, this is because our antenna feed is only sensitive to a limited frequency range, the LNA is only well-characterized for a limited bandwidth, human-made radio interference limits the usable observing band, the ADC only Nyquist samples up to a specific bandwidth before aliasing occurs. All these factors define the usable bandwidth of the instrument. There will be a start frequency $\nu_0$ which is the lowest frequency the system records data at, and a stop frequency $\nu_f$ which is the highest frequency. The range with in those two frequencies is the bandwidth $\Delta \nu = \nu_f - \nu_0$. We can further define two types of bandwidth: analogue bandwidth is the bandwidth of the system based on the analogue filters, and the digital bandwidth which is defines by half the sampling rate of the ADC (the digital bandwidth is set by the Nyquist rate).

The two building blocks of signal filtering is the high-pass filter which allows through signals at high freuqncies while supressing low frequency signals. And, the low-pass filter which is opposite to the high-pass filter. Filter design is big field and a lifetime can be spent on the topic, we will not go into any depth, but just present the basics needed to understand a radio telescope analogue front end.

A filter works by applying a weighted average in the time domain usign some number of samples (called taps). The weights in the average are defined by a window function. A simple window function is a boxcar or truncated Gaussian, but there are many types depending on the desired response. When the time domain signal is transformed to the frequency domain the window function's Fourier transform is called the window response function, the window response of a boxcar is the Sinc function and a Gaussian for a Gaussian.

To set the band of the analogue system we use a bandpass filter which is made up of a low-pass and high-pass filter to created a limted band response. An example of a bandpass filter using a low-pass and high-pass Hann window with 33 taps is shown below. These are finite-impulse response filters (FIRs) which are typically used for bandpass filters, anothe rtype of filter is the infinite-impulse response filter (IIRs) which we will not discuss.


In [ ]:
# Modified from: https://scipy.github.io/old-wiki/pages/Cookbook/FIRFilter.html
# The Nyquist rate of the signal.
sample_rate = 800.0 # ADC sampling rate (MHz)
nyq_rate = sample_rate / 2.0 # Digital Bandwidth

# Low-Pass Filter
# The cutoff frequency of the filter.
lpf_cutoff_hz = 360. # MHz
# Use firwin with a Hann window to create a low-pass FIR filter.
lpf_taps = scipy.signal.firwin(33, lpf_cutoff_hz/nyq_rate, window='hann')

# High-Pass Filter
# The cutoff frequency of the filter.
hpf_cutoff_hz = 40. # MHz
# Use firwin with a Hann window to create a high-pass FIR filter.
hpf_taps = scipy.signal.firwin(33, hpf_cutoff_hz/nyq_rate, window='hann', pass_zero=False)

lpf_w, lpf_h = scipy.signal.freqz(lpf_taps, worN=512)
hpf_w, hpf_h = scipy.signal.freqz(hpf_taps, worN=512)
lpfAmp = np.absolute(lpf_h)
hpfAmp = np.absolute(hpf_h)
bandpassAmp = lpfAmp*hpfAmp

fig, axes = plt.subplots(figsize=(16,6))

plt.subplot(1, 2, 1)
plt.plot((freq0/1e6)+(lpf_w/np.pi)*nyq_rate, bandpassAmp, linewidth=2)
plt.xlabel('Frequency (MHz)')
plt.ylabel('Gain')
plt.title('Bandpass Filter Response')
plt.ylim(-0.05, 1.05)
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot((freq0/1e6)+(lpf_w/np.pi)*nyq_rate, 10.*np.log10(bandpassAmp), linewidth=2)
plt.xlabel('Frequency (MHz)')
plt.ylabel('Gain (dB)')
plt.title('Bandpass Filter Response')
plt.ylim(-10, 0.2)
plt.grid(True)

Design of the bandpass filter depneds on a number of factors, including desired bandwidth, out of band suppression, and cost. Filter design is am engineering challenge set by the project specifications. We will not attempt to cover the topic, see Understanding Digital Signal Processing for an introduction. For our interest, there are a few important characteristics. Usable band: this is the amount of band which can be used for observation, in the example above this is approximately 1175 MHz to 1425 MHz. Roll-off: this is the rate at which the filter drops in response. There is a transtion region in our example from 1100 MHz to 1175 MHz between the bandpass and bandstop regions. The bandpass ripple, we see that the bandpass is not perfectly flat, there is variation, especially when the filter begins to roll-off. We would like a nearly flat response across the band. The filter is continuous in frequency, so even though we would like to put a hard stop to the filter at 1100 MHz and 1500 MHz, the response is non-zero. We will be sampling the analogue signal at some Nyquist rate which will create aliasing effects. We desire to have a low enough filter response at the start and stop frequencies such that signals from out of the band do not alias into the band. The cost of supressing aliasing, is to reduce the overall usable band, or require sharper filters.

The opposite of a bandpass filter is a bandstop filter which supresses signals within a limited frequency range. These are useful in radio astronomy to supress radio frequency interference (RFI) which is strong Human-made radio signals. We discuss RFI in $\S$ 7.8 ➞, for now we need to kow that RFI can be orders of magnitude stronger than astronomical signals. Signals such as mobile telephones at 800 MHz and FM radio aroun 100 MHz (among many others) are present in most populated regions of the Earth.

We can partially seperate out the RFI from the astronomical signals when the observed signal is channelized, but before then RFI can cause issues with the analogue front end due to the dynamic range of the electronic components. Any analogue electornic component is designed to take in a limited range of amplitude as input, this is kown as the dynamic range of the device. Within this dynamic range the device has a linear response, that is the output can be directly mapped to the input based of a linear scaling relation. Input signals out of the dynamic range we produce a non-linear response, and feedback which can affect other components.We would like to build analogue electronics which are optimized to detect astronomical signals, so a solution is to filter out strong RFI signals before these components are used. In the case of strong, persistent, narrow-band RFI (such as FM radio and MObile telephones) a solution is to use a bandstop filter (or notch filter) to supress signals within a limited band around the RFI. Below is a simple spectrum which has three strong, narrow-band RFI source.


In [ ]:
# Spectrum for strong, narrow-band RFI sources
rfiSpectrum = np.ones_like(freqs)
# Add RFI sources
rscale = 20.
idx = np.abs(freqs - 1.33e9).argmin()
rfiSpectrum[idx] += 1. * rscale
idx = np.abs(freqs - 1.34e9).argmin()
rfiSpectrum[idx] += 0.3 * rscale
idx = np.abs(freqs - 1.335e9).argmin()
rfiSpectrum[idx] += 0.5 * rscale

fig, axes = plt.subplots(figsize=(8,6))

plt.plot(freqs/1e9, rfiSpectrum-1)
plt.title('Narrow Band RFI')
plt.xlabel('Frequency (GHz)')
plt.ylabel('Flux (Jy)')
plt.grid(True)

IF we do nothing to supress the RFI the resulting observed spectrum dynamic range will be dominated by the RFI as shown below.


In [ ]:
fig, axes = plt.subplots(figsize=(8,6))

plt.plot(freqs, bandpassAmp * rfiSpectrum, linewidth=2)
plt.xlabel('Frequency (GHz)')
plt.ylabel('Gain')
plt.title('Bandpass with RFI')
plt.grid(True)

Continuing our example, we can design a notch filter to supress 40 MHz of band around these sources. In the figure below the filter is shown in linear and decibel scales. We can see the notch filter will supress signals up to -30 dB in amplitude.


In [ ]:
# Bandstop/Notch Filter
# The cutoff frequency of the filter.
bsf_cutoff_hz = np.array([215., 255.])
# Use firwin with a Hann window to create a low-pass FIR filter.
bsf_taps = scipy.signal.firwin(73, bsf_cutoff_hz/nyq_rate, window='hann')

bsf_w, bsf_h = scipy.signal.freqz(bsf_taps, worN=512)
bsfAmp = np.absolute(bsf_h)

fig, axes = plt.subplots(figsize=(16,6))

plt.subplot(1, 2, 1)
plt.plot((freq0/1e6)+(lpf_w/np.pi)*nyq_rate, bsfAmp, linewidth=2)
plt.xlabel('Frequency (MHz)')
plt.ylabel('Gain')
plt.title('Notch Filter Response')
plt.ylim(-0.05, 1.05)
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot((freq0/1e6)+(lpf_w/np.pi)*nyq_rate, 10.*np.log10(bsfAmp), linewidth=2)
plt.xlabel('Frequency (MHz)')
plt.ylabel('Gain (dB)')
plt.title('Notch Filter Response')
plt.ylim(-35, 1.)
plt.grid(True)

We can then combine the bandstop filter with the bandpass filter, shown below.


In [ ]:
fig, axes = plt.subplots(figsize=(16,6))

plt.subplot(1, 2, 1)
plt.plot((freq0/1e6)+(lpf_w/np.pi)*nyq_rate, bandpassAmp*bsfAmp, linewidth=2)
plt.xlabel('Frequency (MHz)')
plt.ylabel('Gain')
plt.title('Bandpass Filter w/ Notch Filter Response')
plt.ylim(-0.05, 1.05)
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot((freq0/1e6)+(lpf_w/np.pi)*nyq_rate, 10.*np.log10(bandpassAmp*bsfAmp), linewidth=2)
plt.xlabel('Frequency (MHz)')
plt.ylabel('Gain (dB)')
plt.title('Bandpass Filter w/ Notch Filter Response')
plt.ylim(-10, 0.2)
plt.grid(True)

The resulting observed RFI signal (shown below) is significantly supressed due to the notch filter at the cost of losing usable bandwidth. A trade-off that is often necessary.


In [ ]:
fig, axes = plt.subplots(figsize=(8,6))

plt.plot(freqs, bandpassAmp * bsfAmp * rfiSpectrum, linewidth=2)
plt.xlabel('Frequency (GHz)')
plt.ylabel('Gain')
plt.title('Bandpass with Notch Filter for RFI')
plt.grid(True)

Hetrodyne Mixing

After the filters in our block diagram analogue front-end example there is a mixer with a local oscillator. This is know as a hetrodyne detection, which is the process of 'mixing' down the original radio frequencies (RF) $\nu_{\textrm{RF}}$ to an intermediate frequency (IF) $\nu_{\textrm{IF}}$ using a local oscillator (LO) $\nu_{\textrm{LO}}$ such that the signal can be sampled digitally. A local oscillator is a pure, single-frequency tone. Mixing of the RF and LO results in $\nu_{\textrm{IF}} = \nu_{\textrm{RF}} - \nu_{\textrm{LO}}$. Consider two pure sine waves $\sin(2 \pi \nu_{\textrm{RF}})$, $\sin(2 \pi \nu_{\textrm{LO}})$, mixing is the process of multiplying the signals together, using trigonomic identities this results in

$$\sin(2 \pi \nu_{\textrm{RF}}) \cdot \sin(2 \pi \nu_{\textrm{LO}}) = \frac{1}{2} \cos(2 \pi (\nu_{\textrm{RF}} - \nu_{\textrm{LO}})t) + \frac{1}{2} \cos(2 \pi (\nu_{\textrm{RF}} + \nu_{\textrm{LO}})t) $$

Because we have used a bandpass filter, any signal above the stop frequency will be supressed, so the $\cos(2 \pi (\nu_{\textrm{RF}} + \nu_{\textrm{LO}})t)$ term can be dropped, resulting in a mixer output of $\frac{1}{2} \cos(2 \pi (\nu_{\textrm{RF}} - \nu_{\textrm{LO}})t)$. The mixer has shifted the original signal at frequency $\nu_{\textrm{RF}}$ down to $\nu_{\textrm{IF}} = \nu_{\textrm{RF}} - \nu_{\textrm{LO}}$.

The process of hetrodyne mixing is very useful. From an engineering point of view it is much easier to build electronics at lower frequencies than higher frequencies. Using a mixer a high frequency signal is mixed down to lower frequencies, and cheaper electronics can be used. Without used a hetrodyne mixer the alternative is to to use direct RF sampling. This is common is very low frequency radio astronomy (< 200 MHz). As technology improves, direct RF sampling will be used more at higher frequencies, but at the moment most radio telescopes uses hetrodyne mixing.

With nearly all the analogue components in place (we will come to the digitization process soon) we are ready to constructed the observed spectrum of the idealized source (where we started the section) with the various analogue electornics effects included. As we can see in the figure below, it is quite different from the idealized spectrum we started with. The bottom figure shows the time-dependent bandpass change, the right figure is the frequency-dependent gain change.


In [ ]:
# Observed Source Spectrum
# noise spectrum
noiseSpectrum = np.random.normal(loc=0., scale=2., size=gainSpectrum.shape)

# bandpass spectrum
bpSpectrum = np.repeat(bandpassAmp * bsfAmp * rfiSpectrum, ntime).reshape(nchans, ntime)

# observed spectrum
obsSpectrum = bpSpectrum.T * ( (gainSpectrum * np.abs(brightness[:,:,0,0]+brightness[:,:,1,1])) + noiseSpectrum)

# local oscillator frequencies
lo = 1.1e9
ifFreqs = freqs - lo

# save spectrum for use in the next section
# np.savez('data/analogue_spectrum.npz', obsSpectrum, ifFreqs)

fig, axes = plt.subplots(figsize=(8,8))

ax1 = plt.subplot2grid((5,5), (0,0), colspan=4, rowspan=4)
plt.imshow(np.abs(obsSpectrum), aspect='auto')
ax1.get_xaxis().set_visible(False)
ax1.get_yaxis().set_visible(False)

ax2 = plt.subplot2grid((5,5), (4,0), colspan=4, rowspan=1)
plt.plot(ifFreqs/1e9, obsSpectrum[50, :])
plt.plot(ifFreqs/1e9, obsSpectrum[250, :])
plt.plot(ifFreqs/1e9, obsSpectrum[450, :])
plt.xlim(ifFreqs[0]/1e9, ifFreqs[-1]/1e9)
plt.ylabel('Flux')
plt.xlabel('Frequency (GHz)')

ax3 = plt.subplot2grid((5,5), (0,4), colspan=1, rowspan=4)
plt.plot(obsSpectrum[:, 130], np.arange(obsSpectrum.shape[0]))
plt.plot(obsSpectrum[:, 260], np.arange(obsSpectrum.shape[0]))
plt.plot(obsSpectrum[:, 390], np.arange(obsSpectrum.shape[0]))
ax3.get_xaxis().set_visible(False)
ax3.invert_yaxis()
ax3.yaxis.tick_right()
plt.ylabel('Time (s)')
ax3.yaxis.set_label_position('right')

plt.suptitle('Observed Spectrum')

System noise has hidden the smooth source structure. The LNA response has scaled the flux values and added a bit of ripple to the source. The bandpass filter has supressed the source response at the edges of the band. The notch filter has supressed the RFI signal and part of the source signal. The signals range from 0 MHz (DC) to 400 MHz due to the hetrodyne mixer.

Coming back to where we stared the section we can see that the analogue front-end transforms the original astronomical signal into the observed signal. This is what the $\mathbf{G}$-Jones describes. It is generally a time- and frequency-dependent Jones matrix. The time and frequency dependence is due the variation in response from the analogue electronic components used to filter, condition, and amplify the input voltage in order to digitally sample the signal. The generic $\mathbf{G}$-Jones is often broken into a time-stable bandpass $\mathbf{B}$-Jones matrix which represents the filter responses and the general LNA response, and a frequency-stable $\mathbf{G}$-Jones matrix which accounts for the time-dependent variation of the overall analogue system.

7.3.4 Digitization

After signal conditioning in the analogue front-end, we are ready to capture the signal into a digital form so that we can perform the correlation. The term digitization is worth thinking about for a moment, it is the process of taking a continuous-in-time, analogue signal, and converting it to a discrete in time and voltage signal (which are 'digits' described by integers). The last analogue component is the analogue to digital converter (ADC). A simple form of analogue to digital conversion is a cascading resistor ladder which measures discrete steps in an analogue signal, each output of the ladder triggers a bit value of either 0 or 1, thus a digital signal is created. See Chapter 13 of The Art of Electronics.

The dynamic range of an ADC is determined by the number of bits an input analogue signal can be digitized into. For an $n$-bit ADC there are $2^n$ possible digital values, by convention this range is $[-2^{n-1}, 2^{n-1}-1]$. Just as in other electronic components an ADC has input range in which the resulting digitization produces linear results. An astronomical signal needs to be amplified such that the analogue signal is within this linear range. The signal is then digitized, an example of a 4-bit (256 values) ADC is shown in the figure below.


In [ ]:
adcBits = 4 # number of bits in the ADC
adcMin = -1 * 2**(adcBits-1) # minimum ADC value
adcMax = 2**(adcBits-1) - 1 # maximumx ADC value

nsamp = 256 # number of time samples
phi0 = .1*np.pi # starting phase
dphi = 6.*np.pi # length
phis = np.linspace(0, dphi, nsamp)
amp = 6.5 # signal amplitude

analogueSignal = amp*np.sin(phis+phi0)
digitalSignal = analogueSignal.astype(int)

fig, axes = plt.subplots(figsize=(16,8))

plt.subplot(2, 1, 1)
plt.plot(analogueSignal, c='k', label='analogue', linewidth=3, alpha=0.3)
plt.plot(digitalSignal, '.', label='digital', linewidth=3)
plt.ylim(adcMin-1, adcMax+1)
plt.xlim(0, nsamp-1)
plt.axhline(adcMin, ls='--', c='k')
plt.axhline(adcMax, ls='--', c='k')
plt.title('Signal Digitization (4-bit)')
plt.ylabel('Amplitude')
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(analogueSignal-digitalSignal)
plt.xlim(0, nsamp-1)
plt.ylim(-1.2, 1.2)
plt.title('Digitization Error')
plt.ylabel('Amplitude')
plt.xlabel('Time Sample')

The digitization results in discrete-in-time samples, the amplitude of each sample is the closest step in the ADC to the analogue amplitute. This results in digitization error as the conversion to digital results in a loss of information. To reduce the digitization error the maximum number of steps in the ADC should be used. This requires knowledge of the input signal amplitude. If the signal is too weak than only a few bits are used in digitization and the dynamic range is not maximally being used as shown in the figure below.


In [ ]:
adcBits = 4 # number of bits in the ADC
adcMin = -1 * 2**(adcBits-1) # minimum ADC value
adcMax = 2**(adcBits-1) - 1 # maximumx ADC value

nsamp = 256 # number of time samples
phi0 = .1*np.pi # starting phase
dphi = 6.*np.pi # length
phis = np.linspace(0, dphi, nsamp)
amp = 3. # signal amplitude

analogueSignal = amp*np.sin(phis+phi0)
digitalSignal = analogueSignal.astype(int)

fig, axes = plt.subplots(figsize=(16,8))

plt.subplot(2, 1, 1)
plt.plot(analogueSignal, c='k', label='analogue', linewidth=3, alpha=0.3)
plt.plot(digitalSignal, '.', label='digital', linewidth=3)
plt.ylim(adcMin-1, adcMax+1)
plt.xlim(0, nsamp-1)
plt.axhline(adcMin, ls='--', c='k')
plt.axhline(adcMax, ls='--', c='k')
plt.title('Signal Digitization (4-bit)')
plt.ylabel('Amplitude')
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(analogueSignal-digitalSignal)
plt.xlim(0, nsamp-1)
plt.ylim(-1.2, 1.2)
plt.title('Digitization Error')
plt.ylabel('Amplitude')
plt.xlabel('Time Sample')

Another solution to reduce digitization error is to use an ADC with more bits such as the 6-bit ADC figure below. As the number of bits increases the digitization will eventually approach the true analogue signal amplitude.


In [ ]:
adcBits = 6 # number of bits in the ADC
adcMin = -1 * 2**(adcBits-1) # minimum ADC value
adcMax = 2**(adcBits-1) - 1 # maximumx ADC value

nsamp = 256 # number of time samples
phi0 = .1*np.pi # starting phase
dphi = 6.*np.pi # length
phis = np.linspace(0, dphi, nsamp)
amp = 28 # signal amplitude

analogueSignal = amp*np.sin(phis+phi0)
digitalSignal = analogueSignal.astype(int)

fig, axes = plt.subplots(figsize=(16,8))

plt.subplot(2, 1, 1)
plt.plot(analogueSignal, c='k', label='analogue', linewidth=3, alpha=0.3)
plt.plot(digitalSignal, '.', label='digital', linewidth=3)
plt.ylim(adcMin-1, adcMax+1)
plt.xlim(0, nsamp-1)
plt.axhline(adcMin, ls='--', c='k')
plt.axhline(adcMax, ls='--', c='k')
plt.title('Signal Digitization (6-bit)')
plt.ylabel('Amplitude')
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(analogueSignal-digitalSignal)
plt.xlim(0, nsamp-1)
plt.ylim(-1.2, 1.2)
plt.title('Digitization Error')
plt.ylabel('Amplitude')
plt.xlabel('Time Sample')

Now, a problem in radio astronomy is that we wish to maximize the dynamic range of the ADC for a weak astronomical signal but we observe in an environment with strong RFI which will dominate the time-domain telescope signal. We need to sample the full dynamic range of the RFI signal. If the maximum/minimum ADC value is smaller than the RFI signal then the RFI will saturate the ADC, causing the output of the ADC to sit at the highest or lowest value and no sky signal can be captured, see the figure below. This results in spectral leakage where a signal which in reality is localized to a specific frequency 'leaks' into nearby frequencies due to instrumental effects. A strong RFI source is usually a strong sinusoidal wave, when we take the Fourier transform of that signal a stong $\delta$-function is produced in the spectrum and isolated to a single frequency channel. But, if the ADC is saturating then the RFI sinusoid is clipped and turns into a square wave, the Fourier transform of which is similar to the Fourier transform of a boxcar function, a sinc function. A sinc function has longer extended wings which will leak into the other channels.


In [ ]:
adcBits = 4 # number of bits in the ADC
adcMin = -1 * 2**(adcBits-1) # minimum ADC value
adcMax = 2**(adcBits-1) - 1 # maximumx ADC value

amp = 2. # signal amplitude
rfiAmp = 8.5 # RFI signal amplitude
rfiPhi0 = 1.35*np.pi
rfiFreq = 3.37 # frequecy relative to signal frequency

analogueSignal = rfiAmp*np.sin(rfiFreq*(phis+rfiPhi0)) + amp*np.sin(phis+phi0)
digitalSignal = analogueSignal.astype(int).clip(adcMin, adcMax)

fig, axes = plt.subplots(figsize=(16,8))

plt.subplot(2, 1, 1)
plt.plot(amp*np.sin(phis+phi0), c='r', label='signal', linewidth=3, alpha=0.3)
plt.plot(analogueSignal, c='k', label='RFI + signal', linewidth=3, alpha=0.3)
plt.plot(digitalSignal, '.', label='digital', linewidth=3)
plt.ylim(adcMin-4, adcMax+4)
plt.xlim(0, nsamp-1)
plt.axhline(adcMin, ls='--', c='k')
plt.axhline(adcMax, ls='--', c='k')
plt.title('Signal + RFI Digitization  (4-bit)')
plt.ylabel('Amplitude')
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(analogueSignal-digitalSignal)
plt.xlim(0, nsamp-1)
plt.title('Digitization Error')
plt.ylabel('Amplitude')
plt.xlabel('Time Sample')

Not only does saturation cause spectral leakage there is significant digitization error. The naive solution is to attenuate the inout signal such that the RFI does not cause the ADC to saturate. But, if we turn off the RFI we can see in the figure below that the sky signal is not well sampled by the ADC. We are not using the full dynamic range of the ADC. Again, the solution is to use an ADC with more bits to balance RFI from saturating and sampling the sky signal to a sufficient level.


In [ ]:
analogueSignal = amp*np.sin(phis+phi0)
digitalSignal = analogueSignal.astype(int).clip(adcMin, adcMax)

fig, axes = plt.subplots(figsize=(16,8))

plt.subplot(2, 1, 1)
plt.plot(analogueSignal, c='r', label='signal', linewidth=3, alpha=0.3)
plt.plot(digitalSignal, '.', label='digital', linewidth=3)
plt.ylim(adcMin-4, adcMax+4)
plt.xlim(0, nsamp-1)
plt.axhline(adcMin, ls='--', c='k')
plt.axhline(adcMax, ls='--', c='k')
plt.title('Signal Digitization (4-bit, Undersampled)')
plt.ylabel('Amplitude')
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(analogueSignal-digitalSignal)
plt.xlim(0, nsamp-1)
plt.title('Digitization Error')
plt.ylabel('Amplitude')
plt.xlabel('Time Sample')

A typical, modern ADC used for radio astronomy will have 8-12 bits. This is necessary for the strong -RFI environment many telescopes operate in. There is an additional issue for wide-band system, these are telescopes which cover significant fractional bandwidth. Most astronomical radio sources are continuum sources, i.e. the spectrum of the source is smooth and continuous in frequecy. One way to improve the sensitivity of the telescope, as per the radiometer equation, is to increase the observing bandwidth. This is a fine idea, but this means that there is necessarily more power being observed by the telescope, this requires an ADC with increased dynamic range compared to a narrow-band system.

In the past all the electronic components of a radio telescope where analogue, but as we have seen, analogue components add noise and distort signals. Really what we want to do is to turn the observed signal digital. This has the advantage that we can perfectly control the signal at the cost of turning a continuous signal into a discrete one. Analogue components are still required to amplify and filter the signal, but modern radio telescope design aims to convert the signal to digital as soon as possible. Once we have a digital signal we are free to do what we want with the signal, and for interferometry that means correlation.

Future Additions:
  • images: receiver/feed, balun, wave guide, OMT, cryostat, cables, lna
  • sampling rate in ADC section
  • expand radiometer equation discussion
  • expand SEFD discussion, example of observing time for KAT-7, MeerKAT, VLA
  • meerkat or kat-7 analogue receiver chain block diagram
  • full simulation of time-domain voltage signal to auto-correlation spectrum: complex LNA response with phase slope
  • example: KAT-7, MeerKAT waterfall response