Import standard modules:
In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import HTML
HTML('../style/course.css') #apply general CSS
Out[2]:
Import section specific modules:
In [4]:
from IPython.display import HTML
HTML('../style/code_toggle.html')
Out[4]:
TLG:MC: Merge the fundamentals about Fourier series and Fourier transform into chapter two. Keep the physics related to visibilities here, for instance the point source stuff and the Gaussian source stuff. Cut down in size.
TLG:GC: Your English is very good! Well written.
TLG:RC: Reduce your use of subsections. You go way to deep into subsections.
In this section, we will talk briefly about waves, then introduce Fourier analysis. Our aim is to make Fourier analysis as obvious as possible, while still going into sufficient depth to allow students to develop a good physical intuition.
A short note on nomenclature: in the literature, Fourier space can be referred to as phase space, Fourier space, $k$-space, visibility space, uv-space and more besides depending on the discipline. In this book, these terms will be used interchangeably: they all mean the same thing.
TLG:GM: Are these definitions in the glossary? Fourier analysis especially.
TLG:MC: Merge into 2.2. I see what you are trying to accomplish here you want to make a simple intro to a periodic function I would add it to to 2.2 keeping the simplicity though. A periodic function is an endlessly repeating function. The simplest periodic function is the humble sine function:
In [4]:
x=(np.arange(800)-400)/200.*np.pi
plt.title("Sine function")
plt.plot(x,np.sin(x))
Out[4]:
GSF:NC: Add caption and correct numbering
Two quantities define a periodic function: amplitude and frequency. TLG:GM: Are these definitions in the glossary? Amplitude defines the size of the wave; frequency defines how often the wave repeats. In the example above, both are set to 1; the function therefore repeats every $2\pi$ (since the sin function is the y-projection of a particle moving in a circular orbit), and has an amplitude of 1 (since the aforementioned orbit is of unit radius).
Two constants can also be added to a periodic function: amplitude offset and phase offset. TLG:GM: Are these definitions in the glossary? A more general expression for our basic periodic function would therefore be:
In [5]:
amplitude = 1
frequency = 1
phase_offset = np.pi/2.
amplitude_offset = 0
plt.title("General sine function")
plt.plot(x,np.sin(x))
plt.plot(x,amplitude*np.sin(frequency*x+phase_offset)+amplitude_offset)
Out[5]:
GSF:NC: Add caption and correct numbering
I strongly encourage you to play with the graph above - for example, what happens when you set the amplitude to -1 ? What happens if you set the amplitude to -1 and the phase offset to $\pi$?
TLG:ST: Use note styling N.B. Note that setting the phase offset of a sine function to $\frac{\pi}{2}$ gives you a different periodic function: a cosine function. TLG:XX: Use italics for glossary terms only
In many physics textbooks, the frequency $f$ may be referred to as the angular frequency $\omega$. The two are linked by the following relation:
TLG:GM: Are these definitions in the glossary?
TLG:MC: Merge into 2.2. I see what you are trying to accomplish here you want to make a simple intro to a periodic function I would add it to to 2.2 keeping the simplicity though. There exist many periodic functions, all built from the humble sine function. One of particular interest to physicists - for it helps us describe electromagnetic waves - is the exponential form of the complex periodic function, as defined by Euler's formula:
For the remainder of this section, we will use both sine and exponential forms of periodic functions.
TLG:MC: Merge into 2.3. Try to add to that section bringing with it some simplicity. Merge your explanation with what is there. A key discovery of French mathematician Joseph Fourier was that, within a bounded interval, any integrable function - no matter its shape - could be replicated using a linear combination of sines and cosines. The bounded interval is significant: since we are replicating the original function with periodic functions, the predefined interval will essentially correspond to the maximum period of our sines and cosines.
Note that Fourier series are covered in greater mathematical detail in $\S$ 2.3. ➞. Here, we will try to give a more physical understanding of the operation.
Consider some arbitrary function (here, $f(x) = (x-0.5)^3+3$):
In [112]:
x=(np.arange(1200)-600.)/200
y=(x-0.5)**3+2
plt.plot(x,y)
plt.title("Arbitrary function")
plt.ylabel("f(x)")
plt.xlabel("x")
Out[112]:
GSF:NC: Add caption and correct numbering
Let us now consider that we are only interested in this function in the range $x=[-2,2]$. Our interval is thus 4: the period over which we will replicate our initial function $f(x)$ will be 4. Our function here is definitely integrable over the period of interest; we can thus perform a Fourier series and try to replicate our function over our interval of interest. The formula to do so is as follows:
with
Note the factor of 0.5 in front of $a_0$.
In practice, of course, the upper limit of $n\rightarrow\infty$ is not achievable. Let us try to replicate the arbitrary function above using different values of $n_{max}$.
In [120]:
nmax=10
x=(np.arange(1201)-800.)/200
y=(x-0.5)**3+2
def FourierSeriesApprox(xvals,yvals,nmax):
approx=np.zeros_like(yvals)
T=(xvals[-1]-xvals[0])
w=2*np.pi/T
dt=xvals[1]-xvals[0]
approx=approx+1/T*(np.sum(yvals)*dt)
for t in range(len(xvals)):
for n in (np.arange(nmax)+1):
an=2/T*np.sum(np.cos(w*n*xvals)*yvals)*dt
bn=2/T*np.sum(np.sin(w*n*xvals)*yvals)*dt
approx[t]=approx[t]+an*np.cos(w*n*xvals[t])+bn*np.sin(w*n*xvals[t])
return approx
yApprox=FourierSeriesApprox(x,y,nmax)
#plt.plot(x,y)
plt.plot(x,yApprox,label="Fourier series approx. with nmax=%i"%nmax)
plt.plot(x,y,label="Initial function")
plt.legend(loc="upper left")
Out[120]:
GSF:NC: Add caption and correct numbering
Note that our fit becomes worse at the edges if we are fitting a function with a different period than our sampling period. This is because our sines and cosines must repeat at some point - hence the divergence at the edges of our sampling.
Each $n$ corresponds to a given spatial frequency of our function over the sampling range: there are certain frequencies for which our periodic functions prove better "fits" for our function. This is the key point of Fourier analysis: we seek to find what scales contain most of our function's information. As $n_{max}$ increases, we pick up finer and finer detail; as it tends to infinity, we get closer to being able to deal with discontinuities.
Physicists tend to discuss this in terms of energy. Each characteristic scale has a certain amount of information - a certain amount of "energy" from the function's total energy budget. If the function is constant, for example, we know that the characteristic scale will be in the constant term in the Fourier series (i.e. $a_0$): equivalently, we can say that all the function's power is concentrated in $a_0$. The reason we talk about "power" is explained in $\S$ 4.1.2.3 ⤵.
Note that the Fourier series can also be written in terms of a complex exponential as follows:
TLG:MC: Very important but merge into 2.4.
In the example above, we have taken a limited time band and found weights associated with discrete spatial frequencies. However, why should we restrict ourselves to discrete spatial frequencies? Indeed, with sufficiently large $n_{max}$, we should approach a continuous sampling of our function's Fourier coefficients.
By considering our function over its full range (i.e. $[x_0,x_1]\rightarrow[-\infty,\infty]$) and sampling infinitesimally finely in space, we can thus take the complex exponential form of the Fourier series to define the Fourier transform $\tilde{f}$ of a function $f$ defined in Note that Fourier series are covered in greater mathematical detail in $\S$ 2.4. ➞:
In other words, we create a continuous description TLG:GM: Are these definitions in the glossary?, for every frequency, of the relative "weight" of the scale associated with that frequency in our function $f(t)$. For the purposes of this section, that is as far as we will take our analysis; everything that's been said about Fourier series holds for Fourier transforms, which are a more general case of the former.
Note that this works for both time frequencies and spatial frequencies. TLG:GM: Are these definitions in the glossary? This latter quantity is what interests us as astronomers - by taking the (2-dimensional) Fourier transform of a (2-dimensional) image of the sky, we can find where the light of our sources lies. Most importantly, the reverse is true: we can recreate an image if we know its Fourier space well enough! TLG:ST: Use a note style TLG:GC: See the comments about the figures below.
TLG:MC: Very important but merge into 2.4.
TLG:GN: Are F Transform and Series in Glossary? Italic them too.
One key advantage Fourier transforms hold over Fourier series is that we are no longer beholden to poorly-set boundaries. Consider a simple sine wave of the following form:
If we try to find the Fourier coefficient for this function over a range other than a multiple of its period, then we may run into problems. Observe:
In [169]:
def FourierSeriesApprox(xvals,yvals,nmax):
approx=np.zeros_like(yvals)
T=(xvals[-1]-xvals[0])
w=2*np.pi/T
dt=xvals[1]-xvals[0]
approx=approx+1/T*(np.sum(yvals)*dt)
for t in range(len(xvals)):
for n in (np.arange(nmax)+1):
an=2/T*np.sum(np.cos(w*n*xvals)*yvals)*dt
bn=2/T*np.sum(np.sin(w*n*xvals)*yvals)*dt
approx[t]=approx[t]+an*np.cos(w*n*xvals[t])+bn*np.sin(w*n*xvals[t])
return approx
nmax=50
x=(np.arange(1201))/400.
y=np.sin(4*x)
yApprox=FourierSeriesApprox(x,y,nmax)
x1=(np.arange(601-300)*np.pi)/300
y1=np.sin(4*x1)
y1Approx=FourierSeriesApprox(x1,y1,nmax)
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(x,yApprox,label="Fourier series approx. with nmax=%i"%nmax)
plt.plot(x,y,label="Initial function")
plt.legend(loc="upper left")
plt.title("Poorly-set boundary")
#plt.show()
plt.subplot(1,2,2)
plt.plot(x1,y1Approx,label="Fourier series approx. with nmax=%i"%nmax)
plt.plot(x1,y1,label="Initial function")
plt.legend(loc="upper left")
plt.title("Apppropriate boundary")
plt.show()
GSF:NC: Add caption and correct numbering
The only change between those two functions it the $x_{min}$ and $x_{max}$ of our plots - TLG:MC: Rewrite prev sentence. in the first case, we go from 0 to 3, and in the second case from 0 to $\pi$. Although we have only changed our function's boundaries, we get the exact result with $n_{max}=2$ in the second case!. The point here is the characteristic Fourier scales of our function are unchanged - but that boundary errors - poorly-set boundaries, in other words - can cause issues when trying to retrieve these scales.
Fourier transforms do not have this problem: they sample over all space, and therefore never sample under the period of any given characteristic scale of our original function.
Having recovered our function with $n_{max}=2$, and knowing that the boundaries used to do so were $x=[0,\pi/2]$: how many non-zero Fourier coefficients do we have, and what's their value? What do we get with the right boundaries and $n_{max}=1$, and what does that tell us?
We know that the function is odd, so the integral over the range will be zero: $a_0 = 0$. Similarly, because the function is odd, all $b_n$ will be zero, since our function has no even characteristic scales. $a_1$ will be zero, because the function averages to zero over the periods this function investigates. $a_2$, however, will pick up the exact function: it has the correct $w_n$!
Putting a larger $n_{max}$ will change the poorly-bound Fourier series, but not the appropriately-bound one; you can try this for yourself. This means that, in the appropriately-bound case, all $a_{n \ne 2} = 0$ ! Knowing that Fourier transforms never suffer from boundary problems: what is the 1-D Fourier transform of this function? What would it be for $f(x) = \cos(2x) $ ?
As it happens, a sine (or cosine) in the image plane is equivalent to a point in the uv-plane: the Fourier transform of a sine or cosine is a Dirac delta placed at the function's frequency. This is what radio astronomers refer to when they talk about interference fringes: they are talking about the image-plane periodic function associated with a uv-plane point.
TLG:ST: Use a style block
TLG:GC: Here you start moving into the application to inter. I would keep this. TLG:AC: Explain some things here in more detail. Expand section a lot. See figures below. Say assume we have a radio source at the center and assume it is a point source... Mention that there is a link between sources and visibilities. The block diagram. Make the point that there are point sources and how they would look like in visibilities. Do for your examples. The Fourier examples is for the section before this one. You are assuming too much at this point. The students do not know how to link sources with visibilities. Make sure you do not redo proofs in Chap 2. TLG:GN: Point source, Gaussian, dirac delta TLG:AC: Consider addding one baseline sampling the theoretical visibility functions and showing the visibilities. Try to summarize and give some intuition in this first notebook as the rest of the chapter is very theoretical. Consider the extreme case of a point source sitting in the middle of your image. This is typically described by using a Dirac delta - a "pseudo-function" (i.e. a distribution - it behaves, for all intents and purposes, like a function as far as we're concerned here) with the following properties:
Infinitely thin
Infinitely high
Constant area: $\int_{-\infty}^{+\infty} \delta(x)dx = 1$
In other words, you can think of the Dirac delta as the limit where a Gaussian has $\sigma \rightarrow 0$ and $A\rightarrow\infty$, while keeping the area under the Gaussian a constant.
If we have a 1-dimensional image (i.e. a row of pixels) with a point source (i.e. a Dirac delta), what scale is most representative of our image? If we consider its Fourier transform (thus freeing us from questions of boundary), we get the following:
By the definition of the Dirac delta, this immediately gives us:
The Fourier transform of a point source (or, in terms of an astronomical image, of a single-pixel source - an unresolved source) is therefore a constant across all Fourier space. What does this mean in terms of fringes, or spatial frequencies? Quite simply that every fringe TLG:GN: In glossary picks up the point source; it is equally present at all spatial scales. What happens if the Dirac delta is shifted?
Let us consider a point source in a 2-dimensional image. If the source is located at coordinates $(x_0,y_0)$, then the function ${f}(x,y)$ describing the brightness distribution in our image would be:
And the corresponding Fourier space:
Note that here, we define $k_x$ and $k_y$ to be the spatial frequencies in our image, in the $x$ and $y$ directions respectively. Again, by plugging $f(x,y)$ into this equation and using the definition of the Dirac delta, we immediately find:
If our source is at the centre of the image ($x_0=y_0=0$), then we recover our first result: the Fourier plane corresponding to our source is constant. If, however, our source is not in the centre of the image (in radio interferometric parlance, if it is not at phase centre - the reason for this name will be made explicit in later chapters), then our Fourier space consists of the linear combination of two complex waves. In other words, a Dirac delta not located in the centre of the image plane creates fringes in the uv-plane, the periods of which are directly related to its position.
In practice, few sources in the sky are simple point sources. A more "physical" source is a two-dimensional Gaussian:
With different values of $\mu_x,\mu_y,\sigma_x,\sigma_y$, all sorts of "blob-like" sources can be described. More complex sources can be described as a combination of multiple Gaussians. What does the uv-plane of a simple Gaussian source look like? In our case, since $f(x,y) = f(x)f(y)$, we can split the two-dimensional Fourier integral into the product of two one-dimensional Fourier integrals and use the result from Section 2.4.2:
If our Gaussian is located at phase centre (i.e. if $\mu_x=\mu_y=0$), then we are left with another Gaussian:
In other words, the width of our uv-plane Gaussian directly tells us the width of the original Gaussian. If our original Gaussian is located somewhere other than phase centre, we now get a Gaussian multiplied with a fringe pattern in $k$-space (i.e. in the uv-plane). The period of these fringes depends only on $(\mu_x,\mu_y)$ and therefore tell us where the Gaussian is located.
TLG:AC: So you could start this main section listing some Fourier pairs. Point to chap 2 for an expansion.
TLG:AC: ou can use the above block diagram to link an interferometer to the Fourier transform.
TLG:AC: Connect to example in Fourier space
TLG:AC: You might also want to introduce the equation: $V(u,v) = \int I(l,m)e^{-2\pi i (ul+vm)}dldm$.
TLG:MC: Very important material move to a chapter 2.
An important - and often unexplained! - concept in physics is the power spectrum of a function. The name comes from signal analysis. If you receive a given signal - described by some arbitrary function $d(x)$ - then the power distribution of this signal - i.e. literally at what frequencies the power is distributed in this signal - is given by the signal's power spectrum TLG:GN: In glossary:
TLG:GM: Use correct symbol for Fourier Transform. $\mathscr{F}\{\cdot\}$ (also check rest of notebook.
Note that a function's power spectrum is not just its Fourier transform, but rather the Fourier transform of its absolute magnitude squared TLG:GN: In glossary. This means that talking of the "power" placed in a given Fourier coefficient of a function is, strictly speaking, incorrect; however, it is a common and conceptually useful way to talk about the relative values of your Fourier coefficients.
Power spectra are used in a very wide range of applications. Among astrophysicists, they are of particular importance to those studying the Epoch of Reionization.
So far, we have restricted our analyses to images with a single source. Of course, we know that the sky has a lot more than a single source - what then of the Fourier plane of images with more than a single source?
Let us consider an image with $n$ sources of arbitrary shapes, each defined by a functions $f_n(x,y)$. The value of this image at different positions would then be given by:
The corresponding Fourier plane, then, would be the Fourier transform of $B(x,y)$:
In other words, the Fourier plane of our image consists of the linear sum of the Fourier components of each individual source. This means that, in practice, "reading" the Fourier plane of a given image is all but impossible! Some mighty wizards can and do practice such high magic, but they are few and far between.
TLG:ST: Use app style block. Do not use red in the text