The Fourier Transform

 

NOTE to engineers and mathematicians: GET LOST!  This article is a non-rigorous discussion for medical folks.

Above: An aortic pressure signal (minus average pressure) and the sinusoids that were combined to produce the signal.

A frequently employed method in the study of hemodynamics (and many other disciplines as well) is to apply a frequency analysis method called the Fourier transform.  You're likely familiar with the fact that your Doppler ultrasound machine performs an analysis of reflected waves using the Fast Fourier Transform (FFT) to quickly determine (really quickly) the distribution and relative magnitudes of sound frequencies received by the transducer.   "Frequency analysis" of this nature is employed ubiquitously to evaluate signals arising from both natural and artificial sources - obviously in cardiology also.  It's also used for applications such as data filtering and compression.  Just for fun, a jpeg digital image is the result of a two-dimensional Fourier transform that extracts the spatial frequencies of a picture occurring in both the x- and y-axis directions and allows the image to be represented with much less information storage than the original.  This is also to let you know that when talking about Fourier transforms, time is not necessarily the independent variable.

But engineers and scientists employ the Fourier transform in another capacity as well -- to solve mathematical problems.  The use of Fourier methods for studying the cardiovascular system follows naturally due to the  periodic (repetitive) or semi-periodic behavior imposed by the heart rate.  For many purposes it is possible to evaluate the circulation as if it is being subjected to a sinusoidal input; and "sinusoidal" in this context literally means mathematical \(\sin(\omega t)\) and \(\cos(\omega t)\) functions.  If a system is linear, then it turns out that the output resulting from any input is equal to the sum of outputs resulting from any number of individual inputs.  In terms of sinusoids, an input of a given frequency will only cause an output of the same frequency and there will be a constant of proportionality between the two that remains fixed (as long as the system hasn't changed).   A linear approximation like this greatly facilitates analysis and mathematical modeling of the cardiovascular system, even though we know that there are limitations to the approach.

Note that there are a lot of good web sites devoted to this topic and many (most) are more complete and probably more correct than this one.  This one is intended to give you the broad picture, not necessarily doing any of the math yourself, and with a view to cardiovascular-related topics.


Fourier Series

The figure at the top of the article shows how a signal that changes with time, pressure in this case, can be considered as a summation of sinusoidal waves.  It turns out that we could approximate (almost) any function as a sum of sinusoids, but certain types of functions, periodic ones, are particularly suited to it.  Here's what a Fourier series looks like:

\(\Large y(t) = A_0 + \sum_{n=1}^{\infty} A_n \cos\left(\frac{2\pi n t}{T}\right) +  \sum_{n=1}^{\infty} B_n \sin\left(\frac{2\pi n t}{T}\right)   \)

This is a general equation that indicates that we're going to add a bunch of \(\large \cos\) and \(\large \sin\) waves together to approximate (or equal) some function \(\large y(t)\).  It will be the choices for the values of the constants \(\large A_n\) and \(\large B_n\) that determine the actual appearance or character of the function.  \(\large T\) has the same units as \(\large t\) (time?) and constitutes the period of the a repeating function; for the formulation shown you'll get exactly the same value for the function if you stick in some value, \(\large t\), as you would for any other value, \(\large t + n T\), where \(\large n\) is an integer, positive or negative.

The sequence that follows shows how the sinusoids are added in such a way as to reproduce the example pressure tracing.  Note however that sinusoids can be added to reconstruct any signal as long as the signal is Lebesgue integrable and has a finite number of discontinuities (so that an area under the function's curve makes sense), i.e. any biological signal.  The first figure below shows the blue pressure signal as a work in progress as it is reconstructed.  We're simply going to add more and more sinusoids to see how the thing adds up. Initially the blue wave includes only the red average value, and the olive green fundamental frequency.  The latter has a frequency of 1.25 Hz (\(\large T = 0.8\)) corresponding to a heart rate of 75 BPM.  In hemodynamics endeavors, the average value is sometimes referred to as the DC value -- an allusion to "direct current" as in the non-oscillatory current in an electrical circuit.  

For a pseudo-periodic function like this pressure tracing, we only need harmonic frequencies that are integer multiples of the fundamental: 1.25 Hz, 2.5 Hz, 3.75 Hz, 5 Hz,  .. etc., corresponding to n = 1, 2, 3, .. in the above equation.  In the next figure, the 2.5 Hz (purple) harmonic has been added to the blue curve, adding a little bit of detail to its appearance.

Next the 3.75 Hz harmonic ( turquoise) is added. 

As you can see, the addition of each harmonic (5 Hz orange next) adds complexity to the blue pressure signal. We can add as many harmonics as necessary to reproduce the actual signal to the required level of detail.

And finally the original pressure signal is shown with all the harmonics that were necessary to reproduce the pressure tracing. 

 


Some background (If you're math inclined)

Fourier representations of signals are employed ubiquitously in science.  In the following, I'll use \(\large y(t)\) to represent (almost) any function that oscillates in response to independent variable \(\large t\).  Here again is the Fourier representation of a signal \(\large y(t)\):

\(\Large y(t) = A_0 + \sum_{n=1}^{\infty} A_n \cos\left(\frac{2\pi n t}{T}\right) +  \sum_{n=1}^{\infty} B_n \sin\left(\frac{2\pi n t}{T}\right)   \)

Keep in mind, as noted above, that the independent variable doesn't have to be time.  Now it turns out that the sum of any 2 sinusoids of the same frequency is another sinusoid .. of the same frequency.  Also, a sinusoid with a phase shift ( specified by an angle, \(\large \theta\) ), can be expressed as the sum of a \(\large \cos()\) and a \(\large \sin()\) function:  

\(\Large M \cos\left(\frac{2 \pi n t}{T} - \theta\right) = A \cos\left(\frac{2 \pi n t}{T}\right)  + B \sin\left(\frac{2 \pi n t}{T}\right) \)

You can see that the stuff on the right side of the equation is the same as for any one of the frequencies from the Fourier series summation.  The left side represents the same sinusoid, simply expressed differently in terms of a magnitude (or modulus) \(\large M\), and phase shift \(\large \theta\).  So ..

\(\Large y(t) = M_0 + \sum_{n=1}^N M_n \cos\left(\frac{2 \pi n t}{T} - \theta_n\right) \): \(\Large n = 1, 2, 3, ...\)

This says that a periodic function is represented as a sum of sinusoids with period \(\large T\); the specific  function being represented is the result of the values employed for \(\large M_n\) and \(\large \theta_n\).  \(\large M_n\) and \(\large \theta_n\) are the modulus (magnitude) and phase of the sinusoids whose values we have to somehow determine so that \(\large y(t)\) is adequately represented by the sum.   The argument of the cosine function (the bit in the parentheses) ranges over \( \large 2 \pi n \)  (\(\large n\) "revolutions") each time \(\large t\) ranges through a span of the period \(\large T\).   In fact, the pressure signal shown above is reconstructed using a chart of a Fourier series tabulated as modulus and phase, from Hemodynamics by W.R. Milnor.  Here's the table (Note in this example, \(\large -\theta\) has been tabulated and the formula for each sinusoid is \(\large M \cos\left(\frac{2\pi n t}{T} + \theta\right)\):

So the entire complicated pressure tracing above was reconstructed using only a relatively modest amount of information.  The particular waveform requires only a limited number of sinusoids to reconstruct a good representation of the pressure and consists of an average value and 10 sinusoids.  This is a Fourier Series reconstruction where the series contains (only) 11 terms.  

For a given frequency, we need 2 pieces of information to specify the sinusoid; either the magnitude and phase (\(\large M\) and \(\large \theta\)), or separate magnitudes multiplying the \(\large \cos()\) and \(\large \sin()\) functions (\(\large A\) and \(\large B\)). It remains only to show how these quantities are interrelated:

\(\Large A = M \cos(\theta) \)

\(\Large B = M \sin(\theta) \)

and

\(\Large M = \sqrt{A^2 + B^2}\)

\(\Large \theta = \tan^{-1}\left( \frac{B}{A}\right) \)

Note: With \(\large \theta\) defined as show, the following is a trigonometric identity:

\(\Large M \cos\left(\frac{2 \pi n t}{T} - \theta\right) = A \cos\left(\frac{2 \pi n t}{T}\right) + B \sin\left(\frac{2 \pi n t}{T}\right)  = M \cos\left(\theta\right) \cos\left(\frac{2 \pi n t}{T}\right) + M \sin\left(\theta\right) \cos\left(\frac{2 \pi n t}{T}\right)\)  

Some authors choose to represent the Fourier terms as \(\large M \cos\left(\frac{2\pi n t}{T} + \theta\right)\) in which case  \(\large \theta \equiv -\tan^{-1}\left(\frac{B}{A}\right)\).


 To move forward with this topic, we need a fact that was determined in the previous article having to do with how sinusoids and the exponential function are related:

\(\Large e^{j \theta} = \cos(\theta) + j \sin(\theta)\)

It was seen previously that \(\large j = \sqrt{-1}\); this \(\large j\) thing may seem like an idiotic entity, but it turns out to be immensely useful.  In the following, I'll also substitute the angular frequency, \(\large \omega = 2 \pi f\) which has units of radians per second.  Using complex numbers, we can represent an entire sinusoid using only a single complex number \(\large C = A + j\; B\) where \(\large j = \sqrt{-1}\).  When this is done, an expression like ..

\(\Large M \cos(\omega t + \theta) = A \cos(\omega t) + B \sin(\omega t)\)

.. is "simplified" to 

\(\Large C e^{j \omega t} \)

The value of \(\large C\) in this case contains two pieces of information (real and imaginary parts), enough to specify both \(\large A\) and \(\large B\) and therefore  \(\large M\) and \(\large \theta\).  The previous article on complex numbers was all about why this is so.  

\(\Large C e^{j \omega t} = (A + jB) \left[\cos(\omega t) + j \sin(\omega t)\right] = \left[A\cos(\omega t)-B\sin(\omega t)\right] + j\left[B\cos(\omega t) + A\sin(\omega t)\right] \)

As you can see, the real part of \(\large C e^{j \omega t} \) is an entire sinusoid, \( \large \left[A\cos(\omega t)-B\sin(\omega t)\right] \) , that stretches between \(\large -\infty\) and \(\large \infty\). \(\large C e^{j \omega t} \) itself has both real and imaginary parts, but it's taken for granted that we're only interested in the real part in this instance.  (NOTE: All this may seem unnecessarily complicated, but it greatly simplifies working with the math!).  Furthermore, much of the actual mathematical manipulation doesn't even involve the \(\large e^{j \omega t}\) part as suggested by the data table above.  It can be left out until we actually need to produce the sinusoid.   

 Allowing frequency to take on negative values is in keeping with standard pulsed-wave Doppler instrumentation which displays frequencies (velocities) either above or below a 0 baseline, corresponding to towards or away from the transducer.  However the  modulus of the Fourier transform for a real signal is always symmetrical about the y-axis as shown in the above example, so the value at at \(\large +f\) is always identical to the value at \(\large -f\).  (It is an even function).  In contrast, the imaginary part of the spectrum and the phase angle are both odd functions; the spectrum values at \(\large \pm f\) are complex conjugates of each other (see figure below).

The Fourier transform (spectrum) of a real signal (no imaginary component) is even with respect to the modulus and real part (values at \(\large \pm f\) are equal) and odd with respect to the phase and imaginary part (values at \(\large \pm f\) are negatives of each other); the value at \(\large +f\) is the complex conjugate of the value at \(\large -f\).  Complex conjugates are represented graphically above. \(\large a-jb\) is the complex conjugate of \(\large a+jb\) ( \(a\) and \(b\) are real).  The modulus of the complex conjugates are the same; the phases are negatives of each other.


 

The chart above from Milnor / Hemodynamics specifies the modulus and phase for each sinusoid of the pressure signal we've been considering. Although we can think of this as a recipe for construction of the pressure signal from its sinusoids, the chart also enumerates the frequency content - the spectrum of the pressure signal.  Here's the spectrum of the pressure shown as a plot:

 

This gives us a graphical representation of frequency content of the pressure signal and illustrates some typical Fourier (frequency) - domain features of mammalian aortic pressure.  The average pressure value (0 frequency) is the largest Fourier component with a major decrease in the magnitude of subsequent moduli; the oscillatory components are much smaller.  The higher frequencies (e.g. above 6 harmonics) get pretty small and basically cause the little wiggles and bumps in the tracing (the aortic incisura for example).  

I chose a bar graph to suggest that the spectrum contains only specific frequencies - multiples of the fundamental frequency of 1.25 Hz.  In fact, the bars of the graph would need to be infinitely thin to represent the fact that there is absolutely nothing between 1.25 Hz and 2.5 Hz in this signal.   ONLY the discrete frequencies are present.  We'll get more general about Fourier transforms soon, where the signal can have a wide range of frequencies, not limited to specific values.  To put this into familiar terms, we see laminar Doppler signals with only a specific value (a specific Doppler frequency at a particular moment in time) but also spectral broadening where a wide range of frequencies can be present.

Note : There are some good figures and animations on the Wikipedia entry for this topic. 


 Fourier Coefficients

We can represent a function as a sum of sinusoids; all well and good.  But how do we determine the values of \(\large M_n\) and \(\large \theta_n\) ( or \(\large A_n\) and \(\large B_n\)) to use in the series.  These are the coefficients of the Fourier series (Fourier coefficients) which constitute the (frequency) spectrum of the signal. Here's the recipe for determining the Fourier coefficients of a function \(\large y(t)\) with period \(\large T\); the function repeats itself over each time interval of duration \(\large T\):

\(\Large A_0 = \frac{1}{T}  \int_0^T y(t)  dt\)

\(\Large A_n = \frac{2}{T} \int_0^T y(t) \cos\left(\frac{2 \pi n t}{T}\right) dt\)

\(\Large B_n = \frac{2}{T} \int_0^T y(t) \sin\left(\frac{2 \pi n t}{T}\right) dt\)

Oh NO! The dreaded \(\large \int\) thingy!  Hopefully it will take some of the sting out of this to look at some figures so we can get a graphical idea of what this is all about.  First, it's worth pointing out that the \(\large \int\) thing was originally used as a big sloppy \(\large "\text{S}\)" for summation.  We're going to add a bunch of things together.  You're likely aware that an integration is going to come down to an area under a curve.  For the coefficient \(\large A_0\), we multiply the function \(\large y(t)\) by ... \(1\), and determine the area under the curve.  The next figure shows what the math means for the pressure signal we've working on (\(\large y(t)\) is the pressure signal, i.e. \(\large p(t)\)). 

In the top part of the figure, the pressure signal and the function \(1\) are shown.  In the lower part of the figure, the 2 functions have been multiplied together (duh) and the area under the curve is represented in red -- approximately.  I say approximately because you'll see that I've represented the area actually as a lot of skinny rectangles; the height of each approximates the value of the function (\(\large y(t)\)) at each different value of \(\large t\).  I'm going to rewrite the integral now as the approximation shown in the figure:

\(\Large  \sum_{i = 1}^N y_n  \Delta t\)

\(\Large \Delta t = \frac{T}{N}\)

If calculus really bugs you it's worth looking at the last formula closely.  It's the numerical approximation of the expression involving the \(\large \int\) above (the area under the curve) and is most likely the way we'd actually approximate the integral using data digitized from a pressure recording. The \(\large y_n\) are the values of the actual function (\(\large y(t)\)) at the \(\large N\) individual  (discrete) points. The area of each skinny tall rectangle is the product of its base and height, \(\large y_n \Delta t\). The area under the whole curve is then approximated as the sum of \(\large N\) tall skinny rectangles and you can see that the sum of these areas will closely approximate the actual area, depending on how many rectangles are used.  When we allow \(\large N\) to get larger and larger, \(\large \Delta T\) becomes vanishingly small but the approximation become better and better.  In the limit, as \(\large N\rightarrow \infty\) and \(\large \Delta T \rightarrow 0\), the sum becomes the integral.

Now the area under the curve isn't exactly what we want; here's the numerical approximation of the \(\large A_0\) Fourier coefficient.

\(\Large A_0 = \frac{1}{T}  \sum_{i = 1}^N y_n  \Delta t\)

We have to divide the area by the duration i.e. period \(\large T\) which gives us that \(\large A_0\) is the average value (approximately if we're using the summation) of the function over the time interval.  In other words, the value of the Fourier coefficient is the height of a rectangle that has the same area and the same duration as the original function.  It so happens that the average value of the pressure function we've been exploring is approximately \(\large A_0 = 85\) mmHg.  So this was all a very long-winded way of saying that the value of \(\large A_0\) of the Fourier series is just the average value of the function.

But the foregoing will hold us in good stead for the next step.  To find the value of \(\large A_1\) and \(\large B_1\) (Fourier coefficients for the fundamental frequency), we use the same procedure.  The upper part of the next figure shows the pressure function (blue) along with the cosine (red, left) and sine functions (red, right) at the fundamental frequency of 1.25 Hz.  The lower part of the figure shows the function product - what we get by multiplying the \(\large \cos\) (left) or \(\large \sin\) (right) functions by the pressure.  Again the area under the function product is implied by the red skinny rectangles.  Note that area below the 0 baseline is negative and is subtracted from the area summation.

Just as for the average value case, we can sum the area of all the skinny rectangles to determine the approximate area under the curve, and we get some numbers.  Again it's not the area per se that's needed, but the average value so we have to divide by the period \(\large T\).  As you see from the formula, we also have to multiply the terms that have nonzero frequency by 2. After that's all done, it turns out that \(\large A_1 \approx -1.84 \) and \(\large B_1 \approx 18.51 \).  If desired we can then find the value of the modulus, \(\large M_1 = 18.6\) and phase, \(\large \theta = 1.67 \) using the formulas given.

And here's what it looks like to figure out the next Fourier coeffiecient (below).  We just keep stepping through the different frequencies, form the function product (pressure multiplied by the cosine and sine functions), and find its average value over the period of the function. 

Below is the start of the caluculation for \(\large A_8\), showing a relatively high frequency cosine wave (12 Hz) and the function product curve below it.  

 

You can see that the skinny rectangles are getting to the point that they're not all that skinny relative to the cosine wave itself; we have only 10 sample points per cosine wave.  That still turns out to be plenty, but there's a limit to how few times we can sample a sinusoid and not get into trouble.  That limit is the Nyquist limit -- we need a minimum of 2 samples per sinusoid, i.e. we have to sample the signal with at least twice the frequency of the highest frequency sinusoid contained in the signal.  More on that soon.  

Before moving on, and with a view to the next topic, we need to have it in mind that notation for Fourier topics can be variable.  With \(\large T\) being the time interval for the fundamental frequency, \(\large f_0 = 1/T\) is the fundamental frequency  of the signal and \(\large f = n f_0\) represents the range of individual  frequencies in the signal with \(\large n = 1, 2, 3, ..\) taking on integer values.  So we can substitute \(\large f_0\) in place of \(\large 1/T\) in the formulas and \(\large n f_0\) in place of \(\large n/T\).  Lot's of folks (like me) get tired of writing \(\large 2 \pi f\) over and over again and everyone will know what you're talking about if you substitute the angular frequency \(\large \omega\) (introduced above) in place of  \(\large 2 \pi f\).   So we can also substitute \( \omega_0 = 2 \pi f_0\).  These substitutions allow us to write alternate mathematical expressions for the Fourier coefficients and will allow us to progress to a more sophisticated understanding of the Fourier transform.

\(\Large A_n = \frac{2}{T} \int_0^T y(t) \cos(2 \pi n f_0 t) dt =  \frac{2}{T} \int_0^T y(t) \cos(n \omega_0 t) dt  \)

\(\Large B_n = \frac{2}{T} \int_0^T y(t) \sin(2 \pi n f_0) dt =\frac{2}{T} \int_0^T y(t) \sin(n \omega_0 t) dt \)


 Fourier Integral Transform

In the foregoing, it was apparent that we can represent a periodic signal that varies with time as a summation of sinusoids; we used a rather confined set of sinusoids with frequencies that were integer multiples of the fundamental frequency, \(\large f_0 = 1/T\) where \(\large T\) is the period of the signal.  We also saw that the Fourier coefficients, the frequency spectrum of the signal,  is determined from another summation to calculate the area under a curve that we get after multiplying sinusoids by the signal of interest.  We get the spectrum by summing the signal multiplied by sinusoids; we get the signal by summing the spectrum multiplied by sinusoids.  There is a remarkable duality going on here that translates to some even more remarkable mathematical expressions.

In what follows, don't forget our math fact from above:

\(\Large e^{j \theta} = \cos(\theta) + j \sin(\theta)\)

Now in the last section we used only specific sinusoids, ones whose frequencies were integer multiples of a fundamental frequency, \(\large f_0\).  For a Fourier series, we cherry picked the frequencies to use as only the integer multiples of the fundamental frequency.  This time, however, we're going to leave the frequency as a continuous variable in the formulas; the result of the Fourier integral transform will be (potentially) a continuous spectrum, not just the specific values at the predetermined characteristic frequencies.   We'll also allow the frequency  values to range from \(\large -\infty\) to \(\large \infty\).  Don't let the idea of a negative frequency bother you. You've been looking at them for a long time ("positive" and "negative" Doppler shifts depending on whether flow is moving towards or away from the transducer).  While I'm at it, I'll just point out that you can think of the integration from \(\large -\infty\) to \(\large +\infty\) is the "source" of the factor \(2\) in the Fourier series above for the non-zero frequencies.  The integration hits each frequency twice, except  for the the zero frequency. 

Without going into how it was arrived at, here is a definition of the Fourier integral transform:

\(\Large Y(j \omega) = \int_{-\infty}^{\infty} y(t) \; e^{-j \omega t} dt \)  

I'm using lowercase, \(\large y(t)\) to represent the time domain signal and uppercase, \(\large Y(j\omega)\) to represent its Fourier transform; \(\large \omega = 2 \pi f\) as before.  Like the Fourier series, we have an integration of some function (a function of \(\large t\) in this case) multiplied by the \(\large e^{-j\omega t}\) stuff.  But due to the above math fact we have:

\(\Large Y(j \omega) = \int_{-\infty}^{\infty} y(t) \; e^{-j \omega t} dt =  - \int_{-\infty}^{\infty}  y(t) \; \cos (\omega t) \; dt) - j \int_{-\infty}^{\infty} y(t) \; \sin (\omega t) \; dt \)  

Using \(\large e^{-j \omega t}\) in the integrand bundles the \(\large \cos()\) and \(\large \sin()\) parts of the of the Fourier transform into a single compact expression.  The result of the integration will include a real part (multiplying the \(\large \cos()\) function and corresponding to the \(\large A_n\) values that we saw in the above Fourier series) and an imaginary part with the \(\large j\) (multiplying the \(\large \sin()\) function and corresponding to the \(\large B_n\) values above).  Now without further ado, here's the inverse Fourier transform (to sum the frequencies in the spectrum and reproduce the time-domain function):

\(\Large y(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} Y(j \omega) \; e^{+j \omega t} d\omega \)  

I'm always amazed by the similarity and symmetry of the math expressions between the Fourier transform and its inverse.  The first transforms a time-domain signal into a frequency spectrum; the inverse turns the spectrum back into the ("original") time domain signal.   The only differences are the variable of integration ( \(\large \omega\) versus \(\large t\)), the \(\large +\) versus \(\large -\) multiplying \(\large j \omega t\) (and this gets switched back and forth depending on the source),  and the multiplication constant ( \(\large 1/2\pi\)) in front of the inverse transform.  Some sources on this topic don't like the asymmetry and prefer to just fix it, splitting the multiplication constant between the 2 procedures as follows:

\(\Large F(\omega) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} f(t) \; e^{-j \omega t} dt \)  

\(\Large f(t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} F(\omega) \; e^{+j \omega t} d\omega \)  

 So this has been pretty obtuse; let's put it to work and try to get a handle on what it's good for.  We're going to determine the frequency content - the spectrum - of a particular time-domain signal that's relatively easy to do, a step function (other examples are plastered all over the web if you'd care to look):

 

 

This is a signal that is \(\large 0\) everywhere except for the time interval, \(\large -a < t < a\).  You're not going to see anything like this in the circulation, but it so happens that pressure signals like this are used for testing the frequency response of recording gadgets like pressure transducers; we'll see why soon. We can also figure out the frequency spectrum of this signal directly since \(\large y(t)\) in this case is just a constant (let's say \(\large 1.0\)) during the \(\large -a < t < a\) interval and \(0\) everywhere else:

\(\Large Y(j \omega) = \int_{-a}^{a} (1) e^{-j \omega t} dt = \frac{j e^{-j a \omega}}{\omega} - \frac{j e^{+j a \omega }}{\omega}  = j \frac{\cos(- a \omega) + j \sin(- a \omega)}{\omega}  - j \frac{\cos( a \omega) + j \sin(a \omega)}{\omega} = 2 \frac{\sin(a \omega)}{\omega} \)

 The latter expression is the actual frequency spectrum of the time step function. We're not sampling or digitizing the step function -- but using pure math to determine the spectrum exactly.  Here's what the spectrum looks like in a plot with \(\large a = 1\):

This left figure shows both the real and imaginary parts of the spectrum, but we've chosen something of an oddball example where the imaginary part of the spectrum is \(0\) everywhere (due to the symmetry of the step function with respect to the vertical axis).  The modulus of the spectrum is on the right (the phase is either \(0\), when the real part is positive, or \(\pi\), when the real part is negative. You can see that the spectrum oscillates, due to the contribution from \(\large \sin(a \omega)\), but also diminishes with increasing magnitude of the frequency due to having \(\large 2 \pi f = \omega \) in the denominator.  So the step function shown above contains frequencies that range from \(\large -\infty\) to \(\large +\infty\); the magnitude of the contributions diminish with increasing frequency, but they go on forever.  Note: we can't actually compute the value of this function at \(\large \omega = 0\) because of \(0\) in the denominator. However the value approaches a limit: \(\large 2 \frac{\sin(a \omega)}{\omega} \rightarrow 2 a \) as \(\large \omega \rightarrow 0\).

Let's see what happens when we try to reconstruct the time-domain step function from this spectrum using the inverse Fourier transform. In the figure on the left,  10 sinusoids have been added to reconstitute the step function; on the right, 20 sinusoids.

Here's 40 sinusoids (left) and 80 (right).

Here 160!

Here are the points I wanted to make with this:

  1. Yep, the sinusoids described by the spectrum, \(\large Y(j\omega) = 2\frac{\sin(a \omega)}{\omega}\) actually look like they're going to add up to the original step function.  It takes a lot of frequencies to represent a sharp "corner" in the signal; it also takes a lot to represent a long flat region.
  2. Relative to the pressure signal at the top of the article, the step function takes a lot of sinusoids to reconstruct it.  In fact we can see from the formula for the spectrum that it actually goes on forever in frequency.
  3. The reconstructed signal is distorted near a discontinuity.  The step function has 2 discontinuities -- 1) going up and  2) going down.  The reconstructed signal overshoots the actual signal and oscillates for a while before it settles down (known as Gibbs phenomenon).  
  4. This particular example has 2 discontinuities - where the signal goes from 0 to 1.0, or vice versa, over a time interval of 0.  The Fourier reconstruction can't do that.  What it does when enough sinusoids are included, is converge on the value halfway between the 2 values of \(\large y(t)\) on either side of the discontinuity, in this case \(1/2\).  (I don't think there are discontinuities in biological or cardiovascular signals.)

Without proving it, the other thing you should know is that the (whole) Fourier transform of a signal includes ALL the information present in the signal.  i.e. if we have the complete transform, we can get the whole original signal back using the inverse transform procedure.  

Here's the transform for a slightly more general case where \(\large y(t) = 1\) on the interval \(\large a \le t \le b\).  This was just to show we don't always have \(0\) for the imaginary part.

\(\Large  Y(j \omega) = \int_a^b e^{-j \omega t} dt  =  \frac{j e^{-j \omega b}}{\omega} - \frac{j e^{-j \omega a}}{\omega}  =  \frac{\sin(b\omega)-\sin(a\omega)}{\omega}  + j \frac{\cos(b\omega)-\cos(a\omega)}{\omega} \)


 Now we'll turn this thing on it's head and suppose that we have a spectrum - a Fourier transform - that's a step function.   In that case we have  \(\large Y(j\omega) = 1\) on the frequency interval \(\large -a < \omega < a\) and \(0\) everywhere else.  We're saying mathematically that there's some kind of function that has equal quantities (magnitudes) of all frequencies in the low frequency range, but the frequencies just get chopped off at some value, whatever you'd like to specify.   What kind of a time-domain signal has a Fourier transform that looks like that?

\(\Large \int_{-a}^a F(\omega) e^{+j \omega t} d\omega  = \Large  \frac{j e^{-j a t}}{t } - \frac{j e^{+j a t}}{t } = \frac{2 \sin(a t)}{t}\) 

Shoot!  It looks like we get the same dam thing when we determine the Fourier transform of a step function as when we determine the inverse transform of a step function; the only thing different is that one gives a function of time, the other a function of frequency so the \(\large t\) and \(\large \omega\) are interchanged.  This function comes up enough in signal processing that it's been given a name - Sinc.  You've alread seen one, but here's a picture of a normalized Sinc function   \(  \large y(t) = \frac{ \sin(a t) }{a t}\).

 

Normalizing the function does some interesting things.  The function has a peak value of \(1\) at \(\large t=0\).   The area under the curve is also exactly \(1.0\). 


Sinusoidal Functions - Orthogonal and Complete

Here's a little bit about why we can use sinusoids to represent signals.  First, we need to understand the concept of orthogonality of functions.   I'm sure you're used to the idea of orthogonality of vectors.  Two functions, \(\large f(t)\) and \(\large g(t)\) are said to be orthogonal (to each other) on a specified interval, \(\large a \rightarrow b\) if they possess this property:

\(\Large \int_a^b f(t) \; g(t) \; dt = 0\)

The \(\large \sin()\) and \(\large \cos()\) functions possess some very special and interesting properties.  For any given frequency, these 2 functions are always orthogonal to each other over any interval corresponding to an integer multiple of the period.   In the next figure, the \(\large \cos()\) (blue) and \(\large \sin()\) functions are shown (above), and the area under the function product curve is shown below.

 

Hopefully you can imagine that the total red area is zero if the positive and negative areas are added (it is, exactly).   

Not only that; \(\large \cos()\) or \(\large \sin()\)) functions of differing frequency are orthogonal to each other.  The next 2 figures depict \(\large \cos()\) waves of differing frequency:

The lower part of each figure depicts the function product; can you imagine that area components above and below the baseline are identical so the area sum is zero? Over the long haul, the sinusoids are orthogonal even if the frequencies are only slightly different (below):

 

This is one of the necessary properties of sinusoidal functions that allows us to use them to reconstruct signals.  The fact that all the sinusoids are orthogonal to each other guarantees that the Fourier coefficients of a given signal are unique.  

The other necessary property of the sinusoids is that they form a complete set; you can't come up with any signal (that's Lebesque integrable with a finite number of discontinuities) that can't be described by a sum of sinusoids.   It might be an infinite sum, but it will converge on your signal, given enough terms in the series (except at the discontinuities).   You'll have to look online or elsewhere if you want a proof.


Aliasing

Here's something we're all familiar with.  Aliasing comes up whenever we're sampling data at a particular rate that isn't necessarily "fast enough" to represent the frequencies present in the data.  There are lots of misconceptions about this and I'll try not to muddy the waters, but a complete treatment of the topic would require some math that would bore the 2 of you that read this.

Since it was well-said in the Wikipedia entry by someone smarter than me, I'm going to simply quote it:

In general, when a sinusoid of frequency \(\large f\)  is sampled with frequency  \(\large f_s\) the resulting number of cycles per sample is  \(\large f/f_s\) (known as normalized frequency), and the samples are indistinguishable from those of another sinusoid (called an alias) whose normalized frequency differs from  \(\large f/f_s\) by any integer (positive or negative).

In other words (mine this time), the act of sampling data at a constant rate effectively partitions its frequency spectrum into bands (see figure below).  If negative frequencies are excluded,  the bands are \(\large f_s/2\) wide, the width of  the Nyquist frequency.   Sampling a signal of unknown frequency content maps each frequency in the signal to an alias frequency, but makes it impossible to tell which band the frequency arose from.

The Nyquist frequency that you're familiar with, \(\large f_N\),  is twice the sampling frequency, \(\large f_N = 2f_s\).  This is the "highest frequency that can be present in the signal without resulting in alias frequencies".  However, this is an issue that requires some clarification that I'll attempt below.  The alias frequency \(\large f_a\) in terms of the sampling frequency is:

\(\Large f_a = |f - n \;f_s |\)

where

\(\Large n = \text{int} \left[ \frac{f}{f_s} + \frac{1}{2}\right] \)

If we admit the possibility of negative frequencies, as is commonplace for this subject, the formula for the alias frequency is simply:

\(\Large f_a = f - n\;f_s\) where \(\large n\) can be either a positive or negative integer.

The expression for \( \large n\) shows that we're taking the nearest integer portion of \(\large f/f_s\) and \(\large n\) describes the "wrap around" phenomenon that you're familiar with from looking at Doppler signals.  From the appearance of this formula, we obviously would like \(\large n\) to be \(0\) so that the alias frequency and the signal frequency are the same thing.   As soon as the ratio of signal frequency to sampling frequency ( \(\large f/f_s\) ) exceeds \(1/2\), the bracketed quantity exceeds \(1\) and it's the integer portion of the bracketed term that becomes the value of \(n\).  In a nutshell, our problem is that we can no longer determine the value of \(\large n\) after a signal has been sampled.

 

 

 

 

You've probably seen some version of the above figure (hopefully with less nauseating colors).   The color coding shows the frequency bands discussed above (\(\large f/f_s = 0.5\) wide). Sampling maps each frequency in the data to an alias frequency as indicated by the formula and suggested by the figure.  It's said that \(\large f/f_s = 0.5\) is the "folding frequency" because of the way this map appears, layering the actual frequencies back onto the \(\large f/f_s < 0.5\) or \(\large n=0\)  band.    The formula and most texts indicate that all the frequencies are mapped to the \(\large 0 \rightarrow f_s/2\) band - i.e. the purple \(\large n=0\) band in the above figure.  However I'd say that this is simply convention.  Sampling actually destroys any indication of which band the frequencies arose from -- you can no longer learn the value of \(\large n\) depicted in the formulas and figures just given.

Here's a plot of the formula (negative frequencies excluded) -- just another way of showing the figure above.  Whatever frequency is in the signal (x-axis showing  normalized frequency), sampling will map that frequency to the \(\large f/f_s = 0 \rightarrow 0.5\) band (by convention), shown on the y-axis.

The bottom line is that if you wish to preserve the signal subsequent to sampling, you must know something about its frequency content prior to sampling.  To reconstruct the signal, you will need to know enough information about the pre-sampling signal to be able "return" the frequencies to the correct band(s) (the correct value of \(\large n\)).  

Here that is all over again but allowing for the possibility of negative frequencies.  Instead of "folding" the frequencies back onto the \(\large 0 \le f/f_s \le 0.5\) band, the bands are now \(\large 2 f_N\) wide and they're all "stacked" onto the \(\large -0.5 \le f/f_s \le 0.5\) ( \(\large n=0\)) band.

Now some sources will say that you must sample the data at least twice as fast as the highest frequency contained in the signal.  That's a good practice! it will work and it's safe!  Notice that to accomplish this you had to know something about the frequency content of the signal before you sampled it and you took steps to ensure that the signal frequency spectrum is all contained within the \(\large f < f_s/2\) band (\(\large n=0\)). Usual engineering practice is to employ (low-pass) anti-alias filters that ensure negligible frequency content above a cut-off frequency that is less than \(\large f_s/2\).  However this is not a necessary condition. ( Furthermore, we regularly interpret undersampled Doppler signals in cardiology and some examples are shown below.)  If you knew, for example, that all the data arose from the \(\large n=2\) band, you would have the option of sampling it at a much lower rate (called undersampling).  You would be able to determine the actual frequency, \(\large f\), using the formula, the alias frequencies of the sampled data, the sampling frequency, and the value of \(\large n\) that you already know.   You'd simply solve the formula given above for the value of \(\large f\).  That's as far as we need to take that, but there are other circumstances where we'd be able to deduce the actual frequencies from the alias frequencies (the value of \(n\)), given sufficient aprior knowledge of the actual frequency content. 


Now here are some signals to illustrate the concepts. In the first figure below, the (1 Hz) sinusoid is sampled (plot markers) at \(10\) times (10 Hz) the signal frequency; that's enough to describe the sinusoid unambiguously. 

In the next figure, the sinusoid frequency was increased right up to the Nyquist frequency (5 Hz); sampling at 10 Hz, there are exactly 2 samples per wave and this is the minimum required to define the wave unambiguously.

In the next figure, the frequency of the signal (blue, 6 Hz) has increased beyond the Nyquist limit (5 Hz); we could say that the alias frequency appears at 4 Hz as the red sinusoid with sample markers shown to indicate how the 2 waves coincide.  However there an an infinite number of sinusoids that would fulfill the requirement of passing through the sample points; they all just differ by which band they arose from, the value of \(\large n\) in the formula.

 

 In the next 2 figures, the sampling frequency happens to be the same as the signal frequency (10 Hz).  The alias frequency is \(0\) and the alias signal appears (erroneously) as a constant value (red).  Exactly what value depends on the phase relationship between the signal and the sampling.  Both right and left figures are  the same blue signal and sampled (red markers) at the same rate, but with different result.

The final signal below is just grossly undersampled (it's my website; I can put up as many stupid figures up as I want.)

This last figure is here for a reason however.  Although we might say that this high frequency sinusoid was sampled much too slowly to avoid aliasing (it is undersampled), we would still be able to determine the content of the blue sinusoid in the signal as long as we know which band it came from and that there are no other sinusoids in the signal that alias to the frequency of the red sinusoid.

The real problem with aliasing occurs when the frequency content of a signal contains significant quantities (signal power) that alias to the same frequency; at that point in the proceedings we can no longer tell which frequency bands are contributing to the apparent spectrum.   In the figure that follows, the blue signal contains equal magnitude sinusoids at both 1.0 and 11 Hz and sampled at 10 Hz.  Although both of these frequencies have modulus 1.0, the 11 Hz signal aliases at 1.0 Hz and so the sampled data appears to contain only a 1.0 Hz sinusoid with twice the amplitude (the sum of both frequencies).  After sampling, and without knowing the frequency content of the original signal, you're stuck with this.  I don't know of any way to undo this mess unless you knew apriori that the original signal contains equal proportions of these 2 frequencies.  

This is of course what we're stuck with in the case of a pulse wave Doppler recording where there's high velocities (\(\large f > f_N\)) and turbulence.  The mixed frequencies in the turbulence lead to alias frequencies that we have no way of parsing to their appropriate frequency band.

Now that we're used to looking at frequency spectra, let's look at what happens to a signal containing a broad band of frequencies when it gets sampled at a rate that's not fast enough.  As above, a good way to talk about this is to express the frequency as a normalized frequency, dividing the actual frequency \(\large f\) by the sampling frequency, \(\large f_s\).  When we do this, the spectrum is partitioned into bands, each band having a width of \(1.0\).  The "safe" band is between \(-0.5\) and \(0.5\); it contains only frequencies that are being sampled at least twice as fast as the highest frequency in the signal.

For an example, here's a frequency spectrum we've looked at above - the spectrum we'd get for the frequency content of a Step function; the spectrum is described exactly by a Sinc function.  Suppose that the signal is being sampled at a rate suggested by the axis and color coding; everything between \(-0.5\) and \(0.5\) is in the safe zone, everything outside of that zone will be aliased

The next figure indicates where the aliased components of the spectrum will appear.  The frequencies outside the \(\large -0.5 \le f/f_s \le 0.5\) (\(\large n=0\)) band are "stacked" back onto the \(\large n=0\) band so that the apparent spectrum actually appears (erroneously) as indicated by the black line.   

To summarize again:

  1. Once a signal of unknown frequency content has been sampled at a (constant) sampling rate, \(\large f_s\), the entire spectrum of the digitized signal is mapped into a frequency band that is \(\large \pm 0.5 f_s\) in width; if you want to exclude the option of negative frequencies, it's just \(\large +0.5 f_s\) in width.  To say that the frequencies are mapped to the \(\large \pm 0.5 f_s\) band (or the \(\large 0 \rightarrow 0.5 f_s\) band), is arbitrary.  It's simply unknown as to which band(s) contribute to a specific frequency after sampling.
  2. To interpret a spectrum derived after sampling, you must have sufficient knowledge of the original (pre-sampling) signal to interpret the band of origin represented in the spectrum.  The usual statement about this is that you have to sample at least twice as fast as the highest frequency in signal.  That is certainly a sufficient condition, but it isn't necessary.  I would say that you have to know that the spectrum of the original signal is (adequately) confined to a known band.  

 

To put this into everyday context,  next is a PWD aortic flow recording (subcostal, dog) where I've adjusted the pulse repetition frequency (PRF) to a reduced level (reduced the velocity range).  Remember that the vertical scale is really frequency; it's just that the US manufacturers are kind enough to rescale the y-axis for us to read it in terms of velocity.  The two side-by-side images are literally the same recording; they differ only by a baseline adjustment. ( To my eye, peak velocity looks like about 1.3 m/sec)  So, is there aliasing?

And I say, yes there's aliasing, absolutely; and it doesn't matter that I was able to shift the baseline to get the whole flow curve onto the screen.  The figure on the left reveals the true nature of the recording, that the Nyquist frequency corresponds to \(\pm\) ~ 0.8 m/sec.  It shows the signal  wrapping around, in this case returning almost to 0 flow; the latter would correspond to frequency that is twice the Nyquist.   WHY are we able to determine the velocity anyway?

  1. Because you've applied the available filter, your brain.  You know it's an aortic flow, it's all headed away from the transducer, and that the stuff showing up above the 0 baseline is nonsense.  You've chosen which band the frequencies (velocities) are from, even though they beyond the Nyquist limit.  This is an example of undersampling  that you come across regularly.   And ...
  2. The signal isn't contaminated with frequencies from other bands.  We're pretty sure that all of the signal we're seeing is specifically from the band just next to the Nyquist limit (the \( n= \pm 1\) band in terms of the formula for the alias frequency, \(\large f_a = |f - n \; f_s|\). 

 Below is an explanation in terms of the frequency bands that occur when data are sampled.  The left side corresponds to the PW signal on the left above.  We have the velocity range of about \(\large \pm 0.7\) m/sec corresponding to the frequency range \(\large \pm 0.5 f_s\) indicated by the gray band.  However the US manufacturer allows your to apply your "brain filter" and shift the baseline (left PW figure above).  The higher velocities in the signal correspond to a frequency band in the alias range (gray band, \(n=-1\) ).  Although there's aliasing, you know enough about the signal to choose which band the signals arise from. 

We're able to make this determination from Doppler signals (sometimes) because we know that it's flow coming out of a heart and we have some solid expectations as to which way the blood flows, etc.  For the aortic flow, we know that the value starts at zero, just before the valve opens, and progressively increases up to a peak.  So we can actually use the time-history of the signal to help assign the correct frequency band in the presence of aliasing.  The next example is just to make this point, showing an even greater extent of aliasing so that the aortic velocity has wrapped all the way around through 0.  I know you've all seen this at some point in time and I'm NOT recommending the approach.  But we can clearly see on this particular recording that the signal has fully wrapped around once; the peak velocity has been marked on the right at 0.51 m/sec and we can add one full velocity scale, \(\pm\) 1 Nyquist range of 0.8, to get the actual velocity of 1.31 m/sec.  

My point again is that we regularly interpret undersampled signals that exceed the Nyquist limit when we evaluate PW Doppler signals.  This simply requires that we're able to correctly assign the frequency band and that the signal isn't (overly) contaminated with frequencies from more than one band.

 

 

 

RocketTheme Joomla Templates