Thu 18 July 2013 - 15:58:13 CDT
I recently posted a few articles about Fourier decomposition, and I often get asked–and ask myself–“How do I calculate that again?” using a tool like Matlab or (my preference) its open-source cousin, Scilab. Just how we change a data set sampled in time to one sampled in frequency is straightforward. We run the Fast Fourier Transform algorithm:
--> y=fft(x);
Easy. But what if you want to actually understand the data and do something with it? The problems with not knowing how to tailor the output of this calculation are best illustrated with some plots. In Scilab, let’s generate a basic audio signal with two tones sampled at 44.1kHz. First we create an axis consisting of 1024 evenly spaced points in time. A sample rate of 44.1kHz–the CD audio rate–means the sample times are separated by a bit less than 23 microseconds. (Note: this example code works with Scilab 5.4.0 on OSX, your mileage may vary.)
--> N = 1024;
--> fs = 44100;
--> t = (1/fs)*(0:1:(N-1));
Then we create two sine waves with frequencies of 5000Hz and 10000Hz and different amplitudes (1 and 2), add them together, and plot the result.
--> f1 = 10000;
--> a1 = 1;
--> y1 = a1*sin(2*%pi*f1*t);
--> f2 = 5000;
--> a2 = 2;
--> y2 = a2*sin(2*%pi*f2*t);
--> y = y1 + y2;
--> plot(t, y);
After zooming in on the plot, we get something that looks like this:
Taking the Fourier transform and plotting:
--> Y = fft(y);
--> plot(Y);
gives this:
which looks interesting, but it’s not terribly useful. First, we don’t have a frequency axis, just numbered samples. Second, we expect a delta function–a spike (technical term)–for each sine wave, but we see four spikes instead of two plus some negative-going add-ons. And the amplitudes are all messed up; the spikes are too tall.
Second thing first. The result of the FFT is a series of complex numbers, and if you give that to the plot function, Scilab plots the real part by default. In most applications, we’re not interested in real and imaginary parts but rather amplitude and phase. Antennas and image sensors (eyes) and microphones (ears) and speakers and… often only care about intensity, or amplitude-squared; even the phase information is ignored. So for now, let’s plot the amplitude of the FFT using the complex absolute value function.
At the same time, let’s add a proper frequency-axis to the plot. The nature of the FFT is such that the 0th sample is DC or a frequency of 0, and the highest frequency is just less than the sampling frequency. What does “just less” mean? Well, the sample spacing is such that if you extended the sequence by one sample, the Nth point would be the sampling frequency (but because we’re zero referenced, the (N-1)th point is the last one.) An easy way to handle this is:
--> f = linspace(0, fs, N+1); // one sample too many
--> f = f(1:$-1); // drop the last sample
The final trick is to divide by N to scale the amplitude correctly–I always forget that bit, so now I’ve written it down. Now if we plot the scaled amplitude of the FFT vs. our frequency axis
--> plot(f, abs(Y)/N);
we get:
Hmmm, still four spikes instead of two. And we started with amplitudes of 1 and 2, not 0.5 and 1, so the scaling is still wrong. Or is it?
There are two ways to look at this. First, from a Fourier transform table, the transform of a sine function of frequency f is a delta function at f and one at -f each with half the amplitude. So if we move the right half of our plot to the left and shift the axis to include negative frequencies:
--> plot(f-(fs/2), fftshift(abs(Y)/N));
the plot looks like this:
Now we have spikes of amplitude 1–half the original 2–at 5000Hz and -5000Hz, and spikes of amplitude 0.5 at 10000Hz and -10000Hz, just like we expect. (Note that the amplitudes are slightly less than these values, and the spikes flare at their bases. This is a topic for another time, but it’s related to windowing with a rectangle…or convolving with a sinc.)
The other way of looking at the result is to re-examine the unshifted plot and note that the high-frequency spikes are at 39100Hz and 34100Hz. Not coincidentally, these frequencies are the sampling frequency (44100Hz) minus the frequencies of interest (e.g. 44100-5000=39100). The FFT plot is a mirror image folded–origami, get it?–at one-half the sampling frequency. Half the sampling frequency is often called the Nyquist frequency (Harry Nyquist was Swedish, and with folding we’ve resolved this post’s title). The theory says that this is the highest frequency that can be reconstructed from the sampled signal without aliasing. This is also a topic for another time, but here’s a taste.
Let’s keep the 44100Hz sampling rate, but feed it a 30000Hz signal.
--> z = sin(2*%pi*30000*t);
--> plot(f-(fs/2), fftshift(abs(fft(z))/N));
Again we see the characteristic spikes at positive and negative frequencies, but the frequency is not 30000Hz. In fact, it is 44100-30000=14100Hz. So if you record a 30kHz tone at a 44.1kHz rate and play it back, the signal will manifest itself as a tone with a bit less than half the original frequency!
This is fundamental theory that is applied regularly to countless systems: audio, video and imaging, communications, biometric sensors, … While the math and the tools are relatively simple, reality makes the engineering more difficult. Keep on eye on this blog for more.
Mon 20 May 2013 - 15:19:00 CDT
Two previous posts began the discussion of using the Fourier series to break down a periodic signal into constituent sine/cosine waves. Now, to finish that up, we can look at a general case that varies both the phase and the duty cycle of the example square wave. Recall from last time that shifting a square wave by an arbitrary amount would move energy between the constituent cosines and sines, but that we never saw any energy in the even harmonics. The reason for this becomes clear if instead of looking at the Fourier series, we look at the Fourier Transform and employ a couple useful identities.
We use the Fourier Transform to convert a signal between two domains, say time and frequency—although you also often use space and spatial frequency in optics. In this way, it is a generalized form of the Fourier series we have been discussing. The Fourier series breaks down periodic functions into sines and cosines of various frequencies and amplitudes. The Fourier Transform does the same thing for arbitrary signals resulting in amplitudes as a function of frequency.
The two processes also use the same basis although the Fourier transform is often written to use complex exponentials, and the resulting function represents complex amplitudes. But this is just the same thing as the Fourier series! A complex exponential can be written as:
$$ e^{j \theta} = \cos \theta + j \sin \theta $$
…so we can equivalently say that the a coefficients—on the cosine terms—of the Fourier series represent the same thing as the real part of the Fourier transform, and the b coefficients represent the imaginary part.
Back to the square wave. But instead of a square wave, let’s think of this as a combination of three functions. Looking at one period, we have a square pulse, a.k.a. a rect function. A rect function typically ranges from 0 to 1, so the rect is scaled by 2x and it is shifted down. That constant shift is the second function. And finally, the pulse happens periodically, so we can write this as a convolution (lower case) with a Dirac comb, a series of delta functions, infinitesimally thin spikes. Convolving a function with a comb basically makes a copy of the function at each tooth of the comb.
The next figure illustrates this for our square wave. The square wave is in grey; two times the rect (with a width, w, of pi) is in blue; the Dirac comb (with a period, T, of 2 pi) is in red; and the constant offset of -1 is in green.
We can write the function as:
$$ y(t) = 2 \text{rect}\left(\frac{t}{w}\right) \ast \Delta_T(t) - 1 $$
and using some rules of Fourier transforms—and a table of transforms—we can easily write the function’s transform. First, the constant becomes a delta function—one component of a comb—at zero and pointing down because it is negative 1. The rect becomes a sinc function which is a sine divided by its argument, and the argument is given by the width, w, of the rect. The comb becomes a comb with a spacing of 1/T instead of T. For the operators, addition and scaling stay the same, but the convolution becomes multiplication. This—with the rest of the details filled in—gives:
$$ Y(f) = 2 \frac{w \sin(\pi w f)}{\pi w f} \times \frac{\Delta_{1/T}(f)}{T} - \delta(f) $$
These constituents are shown in the top half of the next figure. The sinc (multiplied by 2 and divided by T) is in blue; the comb—not divided by T because we divided the sinc by T for visibility—in red; and the solitary delta function in green. The lower half of the figure shows everything put together. Notice how the resulting function is only non-zero where it coincides with the comb and that the sinc (in grey) represents an amplitude envelope.
Also notice that every other piece of the comb sits at the same frequency as the zero-crossings of the sinc. This explains why we only saw odd harmonics before. The zeros of the sinc sit where the argument of the sine is a multiple of pi or, in this example, with a spacing of 1/w. Because we have a 50% duty-cycle, w is half the period, T, so the zero spacing is 2/T. This is twice the spacing of the comb, so every other harmonic coincides with a zero of the envelope.
Finally, note that the combination of the sinc, comb, and solitary delta at frequency 0 is zero. f=0 represents a constant, or average value—a DC offset if you’re an electrical engineer. Because of the 50% duty cycle, the average value of our square wave is Y(0)=0.
Now let’s wrap up with a picture of what happens when the duty cycle is not 50% and we have an arbitrary shift in the square wave. We’ll have a non-zero average value, and we’ll have a non-zero values for all of the a and b coefficients of the of the Fourier series.
The top half of the plot shows the square wave (black), the sum of the first ten harmonics (blue), and the sums of the first several odd and even harmonics (green and red). The lower half shows the frequency space representation with the amplitudes of the cosine and sine coefficients—or equivalently the real and imaginary parts of the Fourier transform. (The relative amplitudes of the real and imaginary parts are controlled by time-shifting the wave: a time shift is the same as convolving with a single delta function positioned away from t=0, and this is the same as multiplying by a complex exponential—a pure phase term—in the frequency domain.)
So that was kind of a long road over three posts, and it may not seem very practical yet. But understanding transitions between time and frequency and basic Fourier transform properties—especially the convolution-multiplication duality—are fundamental to countless engineering problems in analog and digital signal processing, wireless communication, and optics. (Lenses can do Fourier transforms in space). An upcoming post will look at how understanding these principles helped us find the source of unwanted emissions from a motor drive circuit, so stick around for that.
Thu 09 May 2013 - 11:07:00 CDT
In my earlier post about Mr. Gibbs I discussed using the Fourier series to break down a periodic signal into consitituent sine waves. The example was a 50% duty-cycle square wave, and, as it happens, the Fourier series only included sine terms and only included every other one. Real world signals tend not to be so uniform, and now I want to take a quick closer look at what happens when you break the symmetry.
The first step is straightforward. The earlier square wave was defined as having a transition at t=0 (assuming the signal is a function of time, t). Likewise, sin(t) also transitions through zero at t=0, but cos(t) is has a maximum at 0. So the square wave is in-phase with a sine function, and only exhibits non-zero b coefficients—the coefficients of the sine harmonics—as shown in the following plot.
But if we shift the square wave by one-quarter period, the center of one of its lobes will be at t=0, similar to a cosine. Also, by definition, shifting a sine by a quarter period gives you a cosine. Together this suggests that the shifted square wave is in phase with a cosine, and its Fourier series should have only cosine terms—a coefficients instead of b coefficients. The next plot shows this. And also note that the magnitudes of the coefficients match between the two plots, but due to the mathematical details, the a coefficients alternate between positive and negative values.
Now what happens if the shift is by some amount other than a quarter period (or some multiple)? Remember odd and even functions from high-school math? An odd function follows:
$$ f(-x)=-f(x) $$
A sine function is an example. A cosine function, on the other hand, is even and follows:
$$ f(-x) = f(x) $$
So another way of looking at the two examples above is that the first square wave is an odd function and the second (shifted by a quarter period) is even. If you add even functions together, you end up with an even function. Likewise with adding odd functions to get an odd sum. So it makes sense that the odd square wave only has sine constituents, and the even square wave is built from cosines.
If we do shift the square wave by some arbitrary amount it will no longer be an odd or even function, and it will take some mixture of sines and cosines to build it up. If you shift the sine-phased square wave to the left by some amount t’ and crunch the calculus, you end up with these coefficients:
$$a_n = \frac{2 \sin(n t’)}{n \pi} \left( 1-\cos(n\pi)\right)$$
$$b_n = \frac{2 \cos(n t’)}{n \pi} \left( 1-\cos(n\pi)\right)$$
The next plot shows the result.
When we superimpose the sum of just the cosine terms—the first 10 harmonics—and the sum of just the sine terms on the shifted square wave, a pattern emerges. Note in the plot below that the cosine terms (red) fill in the portion of the shifted square wave that is even, and the sine terms (blue) fill in the portion that is odd.
Now you may have noticed one other detail: regardless of the shift, the coefficients for even harmonics—not talking about even functions, rather simply even indexes (2, 4, 6, …)—are always zero. We’ll look at that in a future post and also look at applying this stuff to some real-world problems. Stay tuned.
Sun 28 April 2013 - 15:55:00 CDT
The math here is driven by what is by far my favorite interview question. In all the interviews that I’ve conducted as an engineer and engineering manager at various places, not one person has ever gotten it right. (So by reading this, you’re giving yourself a leg up.) It’s not a great interview question really, but I just find it interesting that this tidbit stuck with me, and no one else seems to remember it from engineering school.
First, a little background. A great deal of signal processing—and analysis of linear systems in general—relies on basis decomposition. This essentially means taking a signal and breaking it down into a sum of much simpler functions that are each easy to deal with. In Fourier analysis, the basis consists of the set of sinusoids of different frequencies. When adding together the constituent parts, you can vary the amplitude and phase (equivalently the complex amplitude) from one frequency to the next. This is what defines linear analysis/synthesis: basically you can only scale and add.
Further, the individual sines/cosines are called orthogonal—a term borrowed from vector math—because you cannot decompose a sine wave into a sum of sine waves of different frequencies:
$$A \sin(f_1 x) + B \sin(f_2 x) \ne sin(f_3 x)$$
(For unique frequencies and non-zero A and B of course.) Just like if you had three vectors perpendicular (or orthogonal) to each other. In 3D space, no linear combination of a vector in the x-direction and one in the y-direction can result in a vector in the z-direction.
If, in addition, the sine waves forming the basis are scaled appropriately relative to each other, they are said to be normalized. All of this makes an Orthonormal Basis.
Now with that basis—pun intended—if you take a periodic signal, you can find the amplitude of each constituent sine wave making up the signal by calculating the Fourier series. Wikipedia and countless other sources provide the formulas for using sines and cosines to analyze the periodic function f(x):
$$ f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty}a_n \cos(nx) + b_n \sin(nx) $$
where:
$$a_n= \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \cos(nx) dx, \quad n \ge 0 $$
and:
$$b_n= \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \sin(nx) dx, \quad n \ge 1 $$
The coefficients a and b provide the amplitudes of sines (and cosines) with a discrete set of frequencies that are all multiples of the repetition rate of the original periodic signal. These are the harmonics of the original signal. (For non-periodic signals, you calculate the continuous Fourier Transform of the signal to get the amplitudes of a continuum of frequencies of sine waves. Likewise, you can use other basis functions such as the set of Laguerre-Gaussians if you’re dealing in 2D and circular symmetry.)
Again, the good people editing wikipedia have some lovely graphics that you should look at. I will only show a couple basic things here that get us back to the original point. Let’s look at the first three odd, sine harmonics of a square wave. (As it happens, if the square wave crosses zero at the origin and has a 50% duty cycle it only has odd, sine harmonics.) Superimposing these on the square wave, we see that the fundamental (red) has the same period as the square wave itself, and it “fills in” most of the square wave like painting a wall with a big roller. The next two odd harmonics have three (green) and five (blue) times the frequency of the fundamental, and they serve to fill in the corners of the square wave like painting with a trim brush and a detail brush.
Also note that the amplitudes of the harmonics decrease with increasing frequency. Think of it this way: only the sharp transitions benefit from the detail provided by higher and higher frequencies, and these transitions are a relatively small part of the whole signal. And note too that the harmonics and the square wave all cross y=0 together, and the phases of the harmonics alternate in the center of a square pulse (the red and blue traces have peaks where the green trace has a valley).
Now instead of looking at the harmonics individually, we will build the square wave by adding harmonics together. We don’t truly get a square wave until we have added in infinitely many sine waves of increasing frequency. Until infinity—which is never reached—we end up with varying degrees of overshoot and undershoot at the transitions and ripple in the “flats”. The next figure shows the square wave (black) and the sums of the first five (blue), ten (red), and twenty (green) non-zero harmonics.
Now the answer to the interview question: Interestingly, the amplitude of the over/undershoot is constant at roughly 9% of the square wave’s amplitude as you add more and more constituent frequencies. (Rigorously it is not constant, but rather it approaches this finite limit.) This is the Gibbs Phenomenon.
Gibbs Phenomenon. Remember that.
Subscribe to the blog to get automatic update. And post comments or contact us for more information.