13.1 An introduction to the analysis of brain waves  (Page 5/8)

Now we wish to find the coefficients of $f\left(t\right)$ . If we multiply $f\left(t\right)$ by $sin(j·2\pi t)$ and integrate, we get:

Then equating the first and last lines, we have a formula for recovering $a_j$ :

This derivation is valid for as large of $n$ as we want. In fact, the most general application of this is infinite series of sine waves.

Additionally, we can do the same thing with sums of cosines. If $f=b_{1}cos(1·2\pi t)+b_{2}cos(2·2\pi t)+\cdots$ , then we can recover the coefficients in the same way (Exercise 3.1):

And if we have a function that is a mix of the two, such as

we can extract the coefficients by the same equations as above (Exercise 3.2).

Numerical implementation

Now we see why the trapezoid rule was so important. In order to recover the sine coefficients of a function, we need to be able to find the integral of $f$ multiplied against different sine functions. We can do this pretty well with our trapezoid code. To show this method in action, we show the code myfreq.m (full code in Appendix). We'll walk through it here. First, we create a sum of sine waves:

Let's break down the syntax of this operation. We loop through the computation five times, each time representing a frequency j from 1 to 5. Within the loop, the first mytrapz is the integral we care about. We give it two arguments: t and y.*sin(j*2*pi*t) . The first, “ t ", is simply the time grid that we will approximately integrate over. The second corresponds to $f\left(t\right)cos\left(j2\pi t\right)$ , as in equation [link] . Note that j is the frequency we are testing. The operator ( .* ) performs elementwise multiplication, that is, multiplying the corresponding elements in each array. (Many operators such as + and sin() automatically do this.) Thus each entry in the y.*sin(j*2*pi*t) corresponds to one in t .

The second instance of mytrapz() is just normalization (accounts for the length of the wave). After this, we plot the recovered coefficients against the originals. This effectively checks the error in our process.

As we can see from [link] , our method works well.

Frequencies captured by myfreq().

Frequencies captured by matlab's fft().

The remainder of the code uses the Matlab function fft() . This stands for “Fast Fourier Transform." This transforms maps a function from the time domain (the way we normally think about it) to the frequency domain. Basically, we express a signal in terms of the frequency coefficients of which it is composed. The plots from this analysis confirm that our analysis is the same as that of the Matlab-tested functions.