Fourier Stuff - Part 2

Convolution Part 1 Sooooo ....... let's think of how you might go about finding or matching a guess at some signal shape within some sequence of data points that you have. Try this shape, but call it a signal template to be fancy :

yup it's the square wave or box function that we know and love. Here's our signal record :

and you can see just by eye that our template is in there. This would be a fine method for, say, finding the first pulsar ever discovered. Yup indeedy. Jocelyn Bell did it that way on hundreds of yards of paper tape! In our minds we'd probably be doing this :

that is, sliding the template along the record ( or the record past the template ) and saying 'Ah ha!' when they co-incide. But this is way too simple for real work ( though the Stardust@Home project recruited people to do what software probably couldn't ). We need a mathematical, and hence computational/algorithmic technique to (i) slide the template along and (ii) know when a match occurs, or at least assess some degree of match if not exact. That's our motivation. Let's slide the template along but at each offset* note what is the common area under both curves :

here the little gold ball is tracking the leading edge of the moving box. This ball goes up and down according to how much area is simultaneously contained ( green shading ) within both boxes. It hits a peak when they exactly overlap, rises and falls to either side, and there is no common area when the boxes separate completely. Mathematically we would write this as follows :

where I designate x as the offset amount and this is designated as the ( cross ) correlation** ( product ) of the function f with the function g. Note carefully that even though the first named function is included in the product as a complex conjugate ( designated by the trailing asterisk ie. f* ), this will not make a difference for real valued functions. We are accumulating a pointwise arithmetic product for each choice of x. For reasons which will shortly become clearer, I can do what is called a substitution of variable and convert this to :

this operation is called the convolution of f with g. The substitution is a way of adopting a different measurement scale, like feet to metres or gallons to litres, where most importantly the x now appears as a -x. Note that x has units of time in this discussion.

Now I'll do a sidestep ( but not really ). Let's suppose I have the Fourier transforms of f(t) and g(t) respectively, meaning that I have two functions in the frequency domain Ff(s) and Fg(s), and in that domain I'll multiply the functions together as Ff(s)Fg(s) :

this gives a new shape in the frequency domain that is composed of overlapping amounts of each. Note that where either function is zero then the product is zero. I've applied either one of the functions to be the filter of the other one. If I do the (back) transform to the time domain - thus discovering the function that when Fourier transformed gives Ff(s)Fg(s) in the frequency domain - I will obtain :

This is a marvelous result and indeed works both ways. The product of the functions in one domain is the convolution of the corresponding (back) transforms in the other domain. Correlation is equivalent to a convolution, which in turn can be expressed as a multiplication of functions. Most importantly a match/not as judged by correlation in the time domain thus becomes a match/not as judged by the result of multiplication in the frequency domain.

Take an extreme case where f(t) = 0 for all time ie. no signal at all. Ever. Then the transform Ff(s) is zero for all frequencies ( as you'd naturally guess ) because the following integral always returns zero :

... for any frequency choice. [ This is also a case where the Fourier transform gives you the same functional form, like the Gaussian did, but recall that we agreed with eigenvectors to exclude the 'zero eigenvector' instance from eigen-discussions ]. If you multiply this function by Fg(s) which is non-zero in some frequency range ( because g(t) was non-zero for some range in time ) then that gives Ff(s)Fg(s) = 0 regardless of whatever g was up to. This means there is no correlation b/w f and g in the time domain. There is hence no correlation b/w 'no signal' and 'some signal'. While this is intuitively pleasing, there would be no correlation b/w 'no signal' and 'no signal' too ... :-)

Now let's take the other extreme where the template is identically the signal. This should not just be a good match, but in fact the best possible match ( well, there's possibly a normalisation issue ..... ) :

The left hand side of the graphic shows the time domain, the right side the frequency domain. This suggests the general idea a signal match somewhere along the time co-ordinate ought yield some sort of spike within our transform manipulations in the frequency domain. Incidentally for this example we have also 'auto-correlated', giving a measure of total signal presence.

We have now arrived at the point to generally*** describe the core of the E@H 'analysis engine'. Take some stretch of data from Aricebo, Parkes or LIGO IFO's. Do their Fourier transform. You have templates for the signals, and so we do their transforms. Multiply the signal transform by a/the template transform(s) in the frequency domain. Obtain some numerical estimate of match ( ie. to what degree is that function multiple non-zero ) which will hence reflect the degree of match in the time domain. Remember that shifting the signal in time ( slide the template along ..... ) does not affect the magnitude of the Fourier transform at any given frequency.

So why would one bother to do things this round-a-bout way? Why not just convolve ( ie. do the integrals directly in the time domain ) for each choice of template with each stretch of signal data, and at each possible choice of time offset? Efficiency of course is the answer. The Fourier Transform in it's discrete form can be enacted by the FFT**** very quickly, and multiplying two functions together across a single sweep of a common frequency band is not too bad ( and you could integrate/running-total their product as you go ), plus you don't have to transform back from the frequency domain to the time domain to know how good the match is! What one effectively/eventually winds up with goes under the generic heading of 'matched filter' ( matched I think means the filter is specific to the template ), a sort of mathematical gadget that you feed some signal into and out comes a number that reflects the success of the match. It is such numbers that are returned to the E@H servers upon completion of a work unit, in a ranked list of most to least significant 'hits'. This is the 'main crankshaft' of the project.

This makes good intuitive sense. If two time indexed functions are close in shape then they ought have similiar frequency components, and that also implies similiar derivatives. So the degree to which their amounts-per-frequency agree should be close. If you like consider the candidate template as a filter that emphasises the frequencies that match it's own pattern, while diminishing those that don't. An opera tenor singing near a wine glass ( or a grand piano with the lid open .... ) will cause it to vibrate according to the sound energies at the frequencies it's construction will cause it to respond to. The glass will barely vibrate for frequencies off the right pitch. The glass is 'filtering' the tenor's song, and we are mimicking such a physical 'analogue' device but by using numbers within our computers! :-)

Next up : a closer view of convolution in a discrete/integer sense.

Cheers, Mike.

* If I write x' = x + a then f(x) = f(x' - a), so whatever happened at x = 0 now happens at x' = a ( meaning the functional values are the same ). Depending upon your preference this is either shifting the function to the right or shifting the horizontal axis to the left. For x' = x - a then one gets get f(x' + a) and converse shifts.

** Again you may not find all texts agree on nomenclature and exact definitions .....

*** I haven't indicated exactly who is doing what here : server side, client side, pre-processing or post-processing etc. But something akin to the sequence I describe happens somewhere & sometime, no doubt with shortcuts and other cleverness. What also gets rolled into the filter setup/usage is how an interferometer would/should respond to a signal from so-and-so direction in the sky, dopplering due to the Earth's motions, the putative source system's motions, change in frequency at source ...... in addition to what some waveform would 'look like' if just chuffing past you in free space. As I've previously alluded I'm not knowledgeable enough to give further detail on this.

**** We'll see the Discrete Fourier Transform and it's speedy mate the Fast Fourier Transform later on, but here's a quick glance : compare ( back of the envelope ) estimated computational loads via the 'usual' number of 'basic' instructions to execute : a typical DFT transform for N points goes like N^2, but goes like (1/2)N*log_base_2(N) for the FFT. So if one has 'only' 2^10 = 1024 data points, say, then 1024^2 ~ 1 million instructions for DFT but (1/2)*1024*(10) ~ 5000 instructions for the FFT. A factor of around 200 times faster. The relative ratio (DFT divided by FFT) goes like (2N/log_base_2(N)) and hence that benefit increases even more as N increases. If you consider just one second of data taking at a LIGO IFO photodetector that's 16384 = 2^14 points. That translates to a factor of about 8000 times faster with FFT vs DFT!!

So the FFT really does deserve the FAST adjective ..... so much so that one might not even consider doing certain analyses without it .... :-)

Matrix Meandering After some looooong thoughts I've decided to postpone the explanation of integer convolution and discuss ( more ) fully integral/integer indexing of functional values and thus expose the topic of matrices : which we will later see is the 'natural' notation to use when handling the DFT and FFT. Performing a DFT/FFT is really the same as solving a system of simultaneous equations and matrices just rock for that purpose. Somewhat like the complex number business the matrix approach/notation allows simple representation of high level ideas and parks the messy detail under the hood. Plus I reckon it's really cool area of mathematics! The messy detail will turn out to be the inner or dot product, and exactly how one might multiply one matrix by another ( plus whether/why/how/if that operation can be uniquely reverted/inverted ). But first things first. Let's have a look at what restricting stuff to integers as indices looks like on a graph :

so it's a pattern of disconnected dots on a two dimensional plane. The functional values ( vertical distances ) I'm illustrating as integers to be simple, but it's really the fact that the horizontal axis is stepped along in jumps of one unit that matters. We don't actually have a measurement for the in-between bits, and strictly speaking there could be any number and magnitude of functional values within the gaps, for some real world problem. Later we will discuss this very important interpolation issue - or what is the function doing ( probably ) when we aren't measuring it? So this isn't 'smooth' or 'continuous' or 'differentiable' or 'integrable' or any of those concepts discussed earlier. One could write this in tabular form :

which looks like a spreadsheet excerpt. Because it is! That's a good image to have in mind, a matrix is like some part of a spreadsheet. So here's a notation to use to represent these integer situations, called a matrix, starting with a 1-dimensional case or vector :

which one can think of as like a computer language array syntax. A particular number within an array is usually called an array element, with some sort of sub/super-script ( hard to BBCode ). So a_1_2 would denote the element of A in the first row and second column, or a_i_j for the generic i'th row and j'th column element. All the row and column integer indices are positive and traditionally start at one and not zero*. Thus they are the 'counting' or Natural Numbers = {1, 2, 3, 4, 5 .... }. If the values form a single line then the word vector is also applied. I've adopted the convention/assumption that vectors or 1-dimensional matrices have a default vertical column listing of their entries { why I fuss about this detail will become apparent later with matrix multiplication, inner products and ordering of matrices in a matrix product }. But can be written for convenience horizontally by using a little 'T' superscript which stands for ( matrix ) transpose like so :

... or exchange the rows for the columns, a bit like ( ignoring pagination ) English document structure where one goes left to right along a line of glyphs and then go down the lines in order vs. some Asian languages ( I think ) which have consecutive characters going down the columns, upon which finishing a column one moves to the next one on the right ( ? ). Mind you a 1-dimensional matrix can be though of as 2-dimensional : we just have only one row or column, either of which can be any length. As you can see matrices typically have a name or a symbol to know them by and often we will state their dimensions if not immediately apparent. So I can say 'A is an m by n matrix' with the understanding that I mention the number of rows first ( m ) and then the number of columns ( n ). As in : B is a 5 by 3 matrix - so it is called 'B' and has five rows of three columns each. You don't have a different number of columns for each row within a particular matrix. A general matrix is 'rectangular' if m and n are different, and 'square' if m = n ( same number of rows as columns ).

there's a couple of 'obvious' operations with matrices which aren't too surprising :

which works both ways either multiplying in or factoring out a scalar from each and every element. Don't forget we could legitimately use complex numbers, and indeed we will! One can add and subtract matrices but only if each is of exactly the same dimensions, and this is done per each matching row and column entry thus :

..... however matrix multiplication is not similiarly defined ie. multiply the respective elements. Instead it is something that initially looks rather odd. Remember the inner product? For two vectors ( represented as 1-D matrices in this context ) it looks like this :

where I have taken the inner product of a vector with itself to show what that indicates. It is the square of the 'length' of the vector. Good old ( multi-dimensional ) Pythagoras again. The length is also known as the 'norm' or 'magnitude'. We can only sensibly form such a number ( that running total over all index positions ) when the sum makes sense - hence the rule about the number columns on the left equaling the number of rows on the right - which is easily satisfied for a vector and it's own transpose. Next up we will do the biggie ..... an m by n matrix multiplying an n by p matrix ( noting yet again that n is both the number of columns of the first and the number of rows of the second ) which will give an m by p matrix as a result! :-)

Cheers, Mike.

* so the joke goes : hence we are now in the Engineering Department and not the Physics Faculty ..... :-) :-)

Solving An Equation What does it mean to solve an equation? Let us start with what an equation is :

something = another_thing

a statement or declaration that two things are equal. So 5 + 3 = 8 is true but 5 + 3 = 7 is not. All along we have been using the idea of a variable, a named placeholder for some quantity chosen from a certain set or domain. So now I write :

[pre]5 + x = 8 { where x is a real number }[/pre]
and can immediately deduce that x has to be the number 3 for the equation to be true. But if I wrote :

[pre]5 + x = 7 { where x is a real number }[/pre]
then x equals 2 is correct for that. Whereas :

[pre]5 + x = 7.5 { where x is an integer }

or

5 + x = -34 { where x is a natural number }[/pre]
don't have solutions given those choices of domain set. Often we leave out the '{ where x is a real number }' type of qualifying comment as the context has or will define that. But never forget that something like that must be there or is implied to be so. Often a train of logic leading up to an equation being evaluated/solved will have restrictions in place, including the base set from which one draws variable values. Try :

x^2 + 1 = 0

now do we have a value of x which will make this truly an equation ie. left side equals the right? If { where x is a real number } applies then no, but if { where x is a complex number } is applicable then yes. Here is a troublesome 'equation' :

0 = 2

or 'uh oh!' as this is somewhat of a dead end and backtracking is appropriate. No clever choice of domain set will save us here. Often we do reach such cul-de-sac's and have to say 'this is a bad road'. The method of 'proof by contradiction' uses this : assume some road is 'good' -> inevitably reach some statement like 0 = 2 -> thus deduce the road was actually 'bad' -> therefore some other road was the 'good' one ( at least in scenarios where one can enumerate logical outcomes/paths and say 'if this, then not that' ). In any case for '0 = 2' we will say no solution for whatever proposition we were examining that led to that. This is not quite so troublesome :

0 = 0

as in

0 * x = 0

this is not a solution. It is however many the base domain set for x will allow, and that could be an infinite number. This means 'any old x will do' or 'this line of thinking didn't narrow down our possibilities for this variable'. What we did for those uniquely soluble cases above was probably something like this :

5 + x = 8

-5 + 5 + x = -5 + 8

x = 3

that is, we manipulated the equation without doing anything that might break it's equality or truth ( assuming it had those qualities originally ). What we 'actually' did was add these two equations together :

5 + x = 8

-5 = - 5

( I'll often use the word 'equation' without qualification or certainty of truth from now on, but remember the above points of logic ). Thus clearly if :

something_1 = another_thing_1

and also

something_2 = another_thing_2

then

something_1 + something_2

can be evaluated by inserting or substituting the alternate value/expression for 'something_1' and 'something_2'. Hence :

something_1 + something_2 = another_thing_1 + another_thing_2

which I can legitimately do because the equalities pronounced them as interchangeable ie. 'something_1' may swap for 'another_thing_1' and 'something_2' could swap for 'another_thing_2'. So I can in effect add two equations together and not alter their truth*! Similiar logic gives subtraction and also multiplication. Just remember to do to the right side the same operation as suffered by the left**. Division is fine also, but not by zero as that is an invalid operation to perform on either side regardless of their equality or not. This can lead to a subtle trap/error :

[pre]x * y = 10 { x is a real number }

y = 10 / x { x is a non-zero real number }[/pre]
which is generally OK but I have to either know otherwise or make a side note that x = 0 is excluded because of that division. In fact you could have said neither x nor y can be zero from inspecting x * y = 10, by remembering that as 0 is the 'black hole of multiplication' one couldn't get 10 as a product if either x or y was zero. This leads to a more general principle of equations with variables : if one is going to allow the value of the variable to range widely across some set - meaning that any one number in the set can be inserted/substituted into the equation where the variable is mentioned - then it is one's responsibility to check for 'potholes' in the domain set, either existing at the start or introduced during some process. Failing to account for such difficulties/assumptions has been the commonest cause of my math problems anyway. The worst assumption is the one you don't know of ...... :-)

What we aim to achieve is separation. We start with :

something_containing_x = another_thing_possibly_also_containing_x

and get

x = something_not_involving_x

with hopefully the final RHS able to be evaluated without too much trouble like :

x = 4 + ( 7/3 )^4

whatever that turns out to be. Our Fourier Stuff is fortunately all linear equations, indeed I have-been/will-be outlining steps involved in the area of Linear Algebra - an incredibly important topic of maths for E@H and much of science generally. It is when things are non-linear that difficulties can really set in, perhaps intractably. Try this :

cube_root[x*y] = secant_of[PI^y - log(x)]

... best of luck separating that one out for either x or y***! While this equation certainly implies that x and y are related it's not my idea of pleasant mathematical company. :-)

So that's a precis of 'solving an equation'. Next time let's look at a whole gaggle of linear equations. A key point will be how many variables vs how many independent equations do we have in our system.

Cheers, Mike.

* this doesn't apply in Aussie politics .... possibly worldwide politics methinks :-) :-)

** Tread with especial care when taking roots ( and generally inverting functions that are not one-to-one ) as eg.

x^2 + 1 = 0

has two solutions x = i and x = -i

*** I just thought that one up randomly, without any check or analysis, thus I will run the risk of someone actually bothering to attempt to evaluate it .... thus maybe popping out a neat solution! The history of maths is absolutely chock full of guys/gals gnawing on these type of problems for fame or fun or fortune! Or who knows what motivation? If you ever get to read about Donald Coxeter - sadly not with us now - or Andrew Wiles or dear Sir Isaac Newton then you will see what I mean about 'gnawing'. Where would we be without them ? :-)

Well I'm doing what I originally promised : 'as E@H is essentially a part of the LIGO and pulsar signal processing pipelines then I thought I'd explain some of the deeper detail about the mathematical basis of such things ( a sequence of mini-tutes, in fact )'. The audience is whoever wants to read it and the response, if any, is whatever people make of it. Or not. Their choice. I'm not after fame or fortune or testing a textbook chapter I'm writing on anyone. Nothing like that at all. Just for the pleasure of sharing really. And it's a rich vein of knowledge over many centuries that give us what we have today - older guys like me tend to think in such terms. If anything I'm attempting to fill a intellectual void that might be in the mind of some E@H contributors - and those who read the message boards are the more curious minded end of that set - about what those work units are actually doing on their machines. Plus that emphasises why their involvement is such a valuable contribution/gift to the project. Perhaps look upon this as a return to them, other than credits.

But you will note that I do appear to zoom in & out of detail and difficulty level of description, and that is in deference to those who may have done some later high school math but never went on further. So my approach may help to bring them at least some nodding acquaintance with what they would have sampled had they continued. I discovered a while ago that there are some very smart people in the world, way smarter than me ( a guy with four tertiary degrees ), who for better or worse haven't had a 'formal education' for whatever reasons. Maybe I'll meet some more of them by doing this !!

[ More specifically in my experience : dairy farmers. Don't assume that because they look at the back ends of the same bunch of cows twice a day, that they are stupid. The best farmers are terrific integrators of facts. They know what cows are thinking! Try sussing that out with quantum mechanics! Google the gag about 'spherical cows' .... :-) ]

[ I construct these snippets progressively during intervals available between my real/offline life activities. It's a pleasant distraction/contrast. Also I'm enthused, obviously, by this project and seek to diffuse that around. :-) :-) ]

Cheers, Mike.