Intermezzo on Fourier Transforms, Convolution, and Sampling

Fourier Transforms

• Useful reference: R. Bracewell, The Fourier Transform and its Applications. An excellent monograph on this topic. Also includes discussion of the Fast Fourier Transform algorithm, which enabled rapid execution of FTs on the computer.

• The concept of a Fourier Transform (FT) of a function goes back to the idea of trying to represent a function by a sum of sines and or cosines of various frequencies. By having such a sum, one then also knows which frequencies the function contains, so the FT is related to the ``power spectrum'' of a function. The frequencies here are the inverse of the function variable, e.g. a function of time has frequencies in units of s$^{-1}$, a function of position has spatial frequencies in inverse units of length. We talk of ``spatial frequencies'' in that case.

• The significance of FTs is that many FT-relations exist between quantities connected with astronomical observing and data analysis, especially in interferometry, but also in optics of single telescopes, and the theory behind the ``seeing'' for the Earth's atmosphere. In addition, concepts such as FTs, convolutions, (auto)correlations, power spectra, etc play a role in signal and image analysis.

• The general definition of the FT of a function $f(x)$ is:

  
  

  $$F(s) = \int_{-\infty}^{+\infty} f(x) e^{-i2\pi xs} dx$$
  
  
  .

  Note here the link with the cos and sin series, since

  
  

  $$e^{i\theta} = cos\theta + i sin\theta$$
  
  

  The inverse transform is given by:

  
  

  $$f(x) = \int_{-\infty}^{+\infty} F(s) e^{+i2\pi xs} dx$$
  
  
  .

  If you combine these two relations, you express a function f(x) as a double integral of itself. Is this true for all functions? The mathematical conditions are often not met for real signals, but as long as the integral of $\vert f(x)\vert$ over the entire range of $f(x)$ is finite, and any discontinuities in $f(x)$ are finite.

Convolution, Cross-Correlation, Auto-Correlation

• A convolution of one function with another involves a smoothing or averaging of that function, generally resulting in loss of spatial or time resolution. Many processes in obtaining data involve convolutions, for example the blurring of images caused by the atmosphere is a convolution of the original image with a ``point spread function'' that broadens the images. In data analysis, simple actions such as taking running means of data involve a convolution, and the effects of such actions can be easily understood given the connection between convolutions and FTs that we will discuss later on.

• Mathematically, a convolution is defined as follows:

  
  

  $$f(x) \ast g(x) = \int f(u) g(x-u) du$$
  
  

  where the integration is over the entire range in u.

  Convolving thus involves sliding the mirror image of one function along the x-axis, taking the product of that mirrored function and the other function at each position and integrating the area under this product (Draw this process for some simple functions, and visualize what happens!). By replacing g(x) with a simple running mean (e.g. a sequence such as (0.2, 0.2, 0.2, 0.2, 0.2) it is easy to see that this running mean is a convolution of f(x) with this sequence of numbers (in this case replacing each value of f(x) with the mean value of 5 points centered on f(x)). Note: the numbers here are normalized such that their sum is 1 (so the integral $\int_{-\infty}^{+\infty} g(u)du$ is 1. This preserves the overall original scale of values of the function f(x), for example in the case of imaging data, the intensities would be conserved. Since all data is digital, it doesn't hurt to replace ``functions'' by arrays of numbers when thinking about taking convolutions; FTs can also be obtained digitally, through the Fast Fourier Transform algorithm.

• In most cases the convolution functions may be symmetric, in which case it makes no difference whether the function is mirrored or not; however, it is important to recognize that in real life this mirroring does take place in actual situations involving a convolution. This reversal of g(u) to g(x-u) also implies that convolution is commutative; $f(x) * g(x) = g(x) \ast f(x)$. Convolution is also associative, i.e. $f\ast (g \ast h) = (f \ast g)*h$ and distributive over addition, $f*(g+h) = f \ast g + f \ast h$.

• If we do not reverse g(u) to g(x-u), we obtain the cross-correlation of the two functions:

  
  

  $$f \star g = \int f(u) g(u-x) du$$
  
  

  This process is not commutative. Cross-correlating signals allows us to determine how similar the signals are. Likewise, cross-correlation can be used to find sources in data; for example, by cross-correlating the data with an idealized point source function one is in principle able to identify locations and amplitudes of sources as peaks in the cross-correlation map.

• If we take the cross-correlation of a function with itself, we obtain the auto correlation. Like convolution, and cross-correlation, auto correlation involves a loss of information, i.e. it is not possible always to uniquely determine a function if only its auto-correlation is known (the ``phase'' information is lost). Some instruments only measure the auto-correlation of a signal.

• An important process in astronomy involves the attempts to deconvolve signals. This involves trying to obtain higher resolution (``super resolution'') than the images or time series provide, by reversing the convolution process. This only works under certain special circumstances; thinking in terms of FTs, one can only obtain higher frequency information on a function (which is what deconvolution is trying to do) if additional information or constraints on the signal are provided.

• The Energy Spectrum of a function is given by the square of its FT, $\vert F(s)\vert^2$. Just as in the case of an auto-correlation, the spectrum contains no phase information, so there is a loss of information compared to the original function.

Some Special Functions

The relationship between FTs and convolutions is best made clear with the help of some special functions.

• The rectangle function, $II(x)$ is defined such that

  
  

  $$rect (x) = II(x) = 0 for \vert x\vert>1/2, and 1 for \vert x\vert < 1/2$$
  
  

  Question: what would $4II(x/3)$ look like? How about $II(x-8)$?

  This function when multiplied with another function provides an easy way to cut out sections of other functions.

  Question: what would $II(x) \ast II(x)$ look like?

• The sinc function

  
  

  $$sinc (x) = {sin{\pi x}\over {\pi x}}$$
  
  

  This is the filtering or interpolating function; this looks almost, but not quite like a sinusoid of gradually decreasing amplitude. (Plot it!). Its properties include:

  $sinc 0 = 1$

  $sinc n = 0$ where $n$ is a non-zero integer

  $\int_{-\infty}^{+\infty} sinc(x) dx = 1$

  The $sinc$ function is very important. It is the FT of the $II$ function, thus sinc(x) has all frequencies present with equal strength, up to a cutoff frequency. This function, in a convolution, therefore performs perfect low-pass filtering (as we will see when discussing the various theorems). It is also an important interpolation function. If the rectangle function is seen as a single narrow slit (where the width of the slit corresponds to the width of the rectangle function), $sinc^2$ describes the Fraunhofer diffraction pattern of a monochromatic light beam falling through the slit. This illustrates in a simple way the Fourier relation between ``apertures'' and the resulting ``images'' produced by telescopes.

  The two-d analog of the sinc (hence the FT of a square rectangle function in x,y plane) is a Bessel function of the first order. In optics you will learn more about the corresponding analogs of a round 2-d rectangle function (which can be seen as a round telescope mirror).

• Impulse symbol or $\delta$-function. Defined by
  
  $$\delta (x) = 0 for x\not= 0, \int_{-\infty}^{+\infty}\delta (x) dx = 1$$
  
  

  Note: this function is only defined in an integral, and you should not interpret this as $\delta (0) = 1$!

  The role of $\delta(x)$ is quite important; any signal narrower than the resolving power of your instrument or telescope is essentially seen as an impulse function, in a mathematical sense. The function also plays an important role in sampling. The FT of $\delta(x)$ is a function that is 1 everywhere (an infinitely sharp function has all frequencies present at equal strength).

  Question: what is $\int_{-\infty}^{+\infty} \delta(x-a) f(x) dx$?

  Question: what is $\delta(x) \ast f(x)$ ?

• Sampling or replicating symbol, $III(x)$. This is an infinite sequence of unit impulses, or delta functions:

  
  

  $$III(x) = \sum_{n=-\infty}^{+\infty} \delta(x-n),$$
  
  

  and

  
  

  $$III(x) f(x) = \sum_{n=-\infty}^{+\infty} f(n) \delta(x-n)$$
  
  
  .

  Thus, multiplication with $III(x)$ is equivalent to sampling of the function f(x) at regular intervals. How do we control the interval spacing?

  This function is relevant, since in practise we always sample a signal and are never able to measure it completely. Likewise, any reproduction of a singal or data in the computer involves a digitized, sampled representation of the actual data, so all data is sampled. Under convolution, $III(x)$ has a replicating property, reproducing multiple, possibly overlapping, copies of f(x). Try to sketch this yourself. Overlapping will occur if $f(x)$ extends beyond $\vert x\vert>1$. What happens if we convolve $f(x)$ with $III(x/a)$, for different values of $a$?

  The FT of $III(x)$ is $III(s)$ !

Basic Theorems and some other Facts

• We will not prove any of the theorems. If you are curious, Bracewell does this in his book.

• The FT of a Gaussian function is another Gaussian. A narrow Gaussian in one domain will have broad FT, and vice versa. Likewise, any ``narrow'' function will have a ``broad'' FT. Taking this in the limit, one can approximate a the $\delta$ function as an ever narrowing Gaussian of unit integral, and realize that the FT will approach an infinitely broad Gaussian, i.e. a function that is one everywhere, as we stated earlier.

• We can generalize the issue of the relation between the width of one function and that of its FT in the Similarity Theorem:

  If $f(x)$ has FT $F(s)$ then $f(ax)$ has FT ${1\over \vert a\vert} F({s\over a})$.

• Addition Theorem:

  If $f(x)$ and $g(x)$ have FTs $F(s)$ and $G(s)$, respectively, then $f(x) + g(x)$ has FT $F(s) + G(s)$.

• Shift Theorem:

  If $f(x)$ has FT $F(s)$, then $f(x-a)$ has FT $e^{-2\pi ias}F(s)$.

  This implies that a linear shift in one domain corresponds to a phase shift in the FT domain. This makes sense; you are not really changing the frequency content of a function by shifting it, so the amplitude of the FT should not change, only its phase. This theorem finds many examples in optics and radio interferometry; consider a parallel light beam falling on an aperture. To shift the diffracted beam through a small angle, one changes the angle of incidence (i.e. one changes the phase of the illumination across the aperture).

• Convolution Theorem:

  If $f(x)$ has FT $F(s)$ and $g(x)$ has FT $G(s)$, then $f(x) \ast g(x)$ has FT $F(s)G(s)$.

  So, convolution in one domain corresponds to multiplication in the Fourier domain. This theorem makes it often quite easy to visualize the effects of convolutions or FTs, and also to obtain an impression of the FT of complex functions. Often a function can be seen as one or more products or convolutions of simple functions, and successive application of this theorem then makes it easy to visualize what happens in the Fourier domain. Example: an interferometer consisting of multiple dishes can be seen as a single telescope in which many ``holes'' have been cut, or as a replication of many single telescopes. The impact on the subsequent telescope beam (or PSF), which is related to the FT of the aperture plane, can be seen by applying this theorem.

• Properties related to the convolution theorem:
  • The area under a convolution is equal to the product of the areas under the individual functions:

    
    

    $$\int f(x) \ast g(x) dx = [\int f(x) dx][\int g(x) dx].$$
    
    

  • The abscissas of the centers of gravity add under convolution:

    
    

    $$<x>_{f\ast g} = <x>_f + <x>_g$$
    
    

    where $<x>_f = {\int x f(x)dx\over{\int f(x)dx}}$.

  • The second moments add if $<x>_f$ or $<x>_g=0$. In general:

    
    

    $$<x^2>_{f \ast g} = <x^2>_f + <x^2>_g + 2 <x>_f<x>_g$$
    
    
    .

    This has the important consequence that under convolutions, typically the variances add.

• Rayleighs Theorem: The integral of the squared modulus of a function is equal to the integral of the squared modulus of its spectrum:

  
  

  $$\int_{-\infty}^{+\infty}\vert f(x)\vert^2 dx =\int_{-\infty}^{+\infty}\vert F(s)\vert^2ds$$
  
  
  .

  Essentially, this is a statement of energy conservation.

• Autocorrelation Theorem: If $f(x)$ has FT $F(s)$, then its autocorrelation function $\int_{-\infty}^{+\infty} f^\ast(u) f(u+x) du$ has FT, $\vert F(s)\vert^2$. This is a special case of the convolution theorem. So, the autocorrelation of a signal is the FT of its power spectrum. Some instruments work by recording the autocorrelation of a signal rather than the signal itself (in radio astronomy). By taking the FT of that, one has the power spectrum of the signal. Note that phase information is lost in the auto-correlation or power spectrum, so this is not as informative as recording the actual signal.

• Finally, an important property: The definite integral of a function from $-\infty$ to $+\infty$ is equal to the value of its FT transform at the origin:

  
  

  $$\int_{-\infty}^{+\infty} f(x)dx = F(0)$$
  
  
  .

  A consequence of this is that in the case of an interferometer that does not include the ``zero spacing'', hence the 0-meter baseline, the net total recorded flux will be zero. In radio interferometry one speaks of the ``missing short baseline problem.''

Sampling

• Since all digital analysis of signals requires representation of data over finite sections of data (for example pixels on the sky, or time segments in a time series), the question is what the influence of this ``sampling'' of the signal is on the final analysis. In addition, the issue comes up of how frequent the signal needs to be sampled to avoid loss of information, or, if loss is unavoidable, what the nature of that loss of information is.

• Under a certain condition, it is possible to recover the intermediate values (hence the actual signal) with full accuracy. The condition is that the function or signal should be band-limited which means that its FT is nonzero over a finite range of the transform variable, not over an infinite range. In the case of a band-limited signal, the it sampling theorem states the minimum sampling rate required to fully recover the signal:

  A function whose FT is zero for $\vert s\vert> s_c$ is fully specified by values spaced at equal intervals not exceeding ${1/over 2}s_c^{-1}$ save for any harmonic term with zeros at the sampling points.

  Examples: the $sinc$ function is fully specified by the sequence $(x, f(x))$: $....(-2,0), (-1,0), (0,1), (1,0), (2,0), ....$ even though all except one sample are 0! On the other hand, the sequence $....(-2,0), (-1,0), (0,0), (1,0), (2,0),...$ does not fully specify the function $sin(\pi x)$, because of the exception regarding harmonic terms stated above.

• In practise, a function or signal may not strictly be bandlimited. However, there is usually a range in $s$ beyond which $F(s)$ becomes negligibly small; in that case we can characterize the error or loss of information resulting from the imperfect sampling. In this case, the measuring device may impose the cutoff frequency, so that the FT becomes $II({s\over{2s_c}})G(s)$. Here $G(s)$ is the FT of the actual signal and $s_c$ is the cutoff frequency imposed by the measuring device (e.g. the pixel size on a CCD restricts the rate of sampling of the point spread function in imaging).

• Since the FT of $II({s\over{2s_c}})G(s)$ is $sinc2s_cx \ast g(x)$, the effect imposed by the sampling is to convolve the original signal with a sinc function.

• How we recover the values of the original signal at intermediate, non-sampled points is the subject of interpolation. Interpolation is extremely common in astronomical data analysis. For example, one often has to re-grid an image to align it with another image. This re-gridding could involve a change in pixel size, image rotation, and image translation. It is important to recognize that interpolation is rarely perfect, i.e. usually it implies further smoothing (that is, loss of resolution) of the sampled signal. For a band-limited signal or function, convolution with $2s_csinc(2s_cx)$ will yield values for all x which are exact if the sampling theorem was met.

• If the function is not band-limited, (e.g. Gaussian) then the $sinc$ function is not necessarily the best interpolation function; a Gaussian would be better. See .e.g Press etal's monograph on Numerical Recipes for discussion of various interpolation methods.

• In practise, any interpolation will be done by using only a finite range of points in x, so even using a $sinc$ interpolation function in practise will not produce the ideal result unless you take this $sinc$ along over a very wide range in x. It is of course possible to estimate the consequences of this. Only including the $sinc$ over a limited range in $x$ implies that you are multiplying the $sinc$ with a rectangle function, so using FT theorems, one may visualize what the consequences can be.

• Undersampling. What happens if you undersample a signal or function? See the figure I handed out. (Figures 10.3 and 10.6 in Bracewell). Essentially, undersampling a function causes overlap in the frequency domain, so the higher frequencies are messed up; this is called aliasing. Figure 10.6 demonstrates the effect of sampling a function at too low a rate to detect a spike. In astronomy, the situation in e.g. imaging is slight different, since the pixels have finite extent, and adjoin smoothly. The consequence of this is that even if we do undersample (by having pixels too large for the point spread function), we do not fully miss stars since the pixels adjoin; there is flux conservation, but what happens is that the stellar profiles are broader than they could be, and stars may overlap. We are still undersampling the signal though.