# Connexions

You are here: Home » Content » The Art of the PFUG » An Introduction to the Analysis of Brain Waves

• How to Compose a PFUG

### Lenses

What is a lens?

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

#### Affiliated with (What does "Affiliated with" mean?)

This content is either by members of the organizations listed or about topics related to the organizations listed. Click each link to see a list of all content affiliated with the organization.
• Rice Digital Scholarship

This collection is included in aLens by: Digital Scholarship at Rice University

Click the "Rice Digital Scholarship" link to see all content affiliated with them.

#### Also in these lenses

• Lens for Engineering

This collection is included inLens: Lens for Engineering
By: Sidney Burrus

Click the "Lens for Engineering" link to see all content selected in this lens.

### Recently Viewed

This feature requires Javascript to be enabled.

Inside Collection (Book):

Book by: Steven J. Cox. E-mail the author

# An Introduction to the Analysis of Brain Waves

Module by: Ryan George, Steven J. Cox. E-mail the authorsEdited By: Ryan George, Steven J. Cox

Summary: This module builds the tools necessary for the frequency analysis of brain waves as recording by an electroencephalograph. We proceed from the Pythagorean Theorem to sine waves, the trapezoid rule and finally to Fourier decomposition.

0. Background

1. Sine and Cosine Waves

2. Trapezoid Rule for estimating area

3. Fourier Method for decomposing signals

4. Spectrogram application to analyzing brain waves

## Background: Brain Waves and the EEG

Signals are sent through the brain using both chemical and electrical means. The synchronized electrical activity of individual neurons adds up to something big enough to detect on from outside the head. To measure it, we use a set of electrical nodes called an electroencephalogam (EEG). The measured activity reflects different states of the brain which in turn tell us something about the mindset of the person. Our goal in this module is to decompose an EEG signal into its different frequencies, which is intuitively the most meaningful piece of information.

## Sine and Cosine Waves

Brainwaves have complex shapes that are not easily interpreted. In order to study these waves, we need to develop some mathematical tools that will tell us about different waves. To outline, we begin by talking about pure (sine or cosine) waves, then move to the trapezoid rule for estimating area under a curve. Next, we develop Fourier analysis for picking out the frequencies in a jumbled signal, and finally use these tools to create spectrograms, which allow us to track different frequencies over time.

### Sine Waves

The sine wave is a mathematical function. It describes many physical phenomena, including sound waves and oscillation. It looks just like a wave. MATLAB uses the sin function to make sin waves. For example, to make Figure 1, we use the code:

>>t = 0:.01:1;
>>y = sin(2*pi*t);
>>plot(t,y);


The sine wave is defined by the lengths and angles of a triangle. Run sincirc.m (copied below) to see how the sine and cosine values relate to the angle ϕϕ of the triangle. As you can see, if ϕϕ is the angle of a right triangle with hypotenuse 1 (illustrated by the circle) , sin(ϕ)sin(ϕ) is the height of the triangle and cos(ϕ)cos(ϕ) is the base of it:

% sincirc.m
%
% sincirc.m illustrates the relation of the sin and cosine waves to the circle.

%define parameters
Nturns = 2;
steps_per_turn = 9;
step_inc = 2*pi/steps_per_turn;

%set up points for circle
circ_x = cos(0:.01:2*pi);
circ_y = sin(0:.01:2*pi);
axis equal

%loop over triangles with different angles
for n = 1:Nturns * steps_per_turn;
phi = n * step_inc + pi/4;

%plot circle, then triangle, then text
plot(circ_x, circ_y);
axis([-1 1 -1 1] * 1.5);
line([0 cos(phi)], [0 sin(phi)]);
line([1 1] * cos(phi), [0 sin(phi)]);
line([0 cos(phi)], [0 0]);
text(cos(phi)/2 , -.1*sign(sin(phi)),'cos(\varphi)')
text(cos(phi) + .1*(sign(cos(phi))-.5), sin(phi)/2, 'sin(\varphi)')
text(cos(phi)*.2, sin(phi)*.1,'\varphi');
pause(.5);
end


### Characteristics of the Sine Wave

The sin wave has three primary characteristics:

#### 1.

Frequency measures how often a wave passes. We can make a wave with frequency ωω by writing:

sin 2 π ω t sin 2 π ω t
(1)

Aside: The wave described by sin(t)sin(t) has frequency 1/2π1/2π. If we instead write sin(2πt)sin(2πt), we will have a wave with frequency 1, which is easier to work with.

We can express the same information in terms of wavelength. Wavelength is how close neighboring waves are to each other. It is inversely proportional to frequency, which means that the higher the frequency, the smaller the wavelength. If is wavelength, we write this in the following way:

= 1 ω = 1 ω
(2)

A wave with wavelength is written as sin(2πt)sin(2πt).

#### 2.

Amplitude measures how high the wave is. We can make a wave of amplitude aa by writing:

a · sin ( 2 π t ) a · sin ( 2 π t )
(3)

Sometimes aa will be negative. In this case, we still say that the wave has amplitude aa, but note that the function will be flipped across the x-axis.

#### 3.

Phase describes how far the wave has been shifted from center. To create a wave with phase pp, we write:

sin ( 2 π ω t + p ) sin ( 2 π ω t + p )
(4)

Since a sin wave repeats every 2π2π, the following is always true (that is, for any ϕϕ):

sin ( ϕ ) = sin ( ϕ + 2 π ) sin ( ϕ ) = sin ( ϕ + 2 π )
(5)

Aside: A cosine wave is a sine wave shifted back by a π/2π/2, a quarter of the standard wave:

cos ( ϕ ) = sin ( ϕ + π / 2 ) cos ( ϕ ) = sin ( ϕ + π / 2 )
(6)

Every sine (cosine) wave can be described completely by these characteristics. These are shown in Figure 2 (phase =0=0 for simplicity):

To code Figure 2 in MATLAB, use:

>>t = 0:.01:1;
>>amp = 3;
>>freq = 1;
>>phase = 0
>>y = amp*sin(freq*2*pi*t + phase);


#### Hearing Sine Waves

The sine wave represents a pure tone. To hear one, we use the MATLAB function sound(), which converts a vector into sound. Find the frequency that is specified, and compare to human range of hearing. Should we be able to hear this sound?

>>freq = 1000;
>>samp_rate = 1e4;
>>duration = 1;
>>samples = 0 : (1/samp_rate) : duration;
>>sound_wave = sin(2 * pi * samples * freq);
>>sound(sound_wave, samp_rate);


Enter the above code into Matlab to hear the sound.

We can add together multiple sin waves to accomplish different shapes. For example, if we add the wave [4·sin(2πt)][4·sin(2πt)] and the wave [sin(6·(2πt))][sin(6·(2πt))] we get:

Look at the figure and try to identify the effect of each wave. The first wave has frequency 1 and amplitude 4. This accounts for the large up-and-down motion that only goes through one cycle in the figure. The second wave has frequency 6 (wavelength 1616) and amplitude 1. This accounts for the small wiggles that happen many times in the figure.

Figure 4 is implemented in MATLAB with the following code:

>t = 0:.01:1;
>>y = 4*sin(2*pi*t)+sin(6*2*pi*t);
>>plot(t,y)


A wave of any shape can be expressed as a sum of sin and cosine waves, although it may take infinitely many. In the end of this module we will find interesting uses for the this fact.

### Exercises

#### 1.1

Figura 4 two shows a sum of two sine waves (with phase 0). Try to replicate it, and report what the amplitudes and frequencies of each wave is.

#### 1.2

Use the identities cos(ϕ)=sin(ϕ+π/2)cos(ϕ)=sin(ϕ+π/2) and sin(ϕ)=sin(ϕ+2π)sin(ϕ)=sin(ϕ+2π) to solve for pp in the following equations by making the above substitutions. Check your answer visually against the figures of sine and cosine waves:

#### Equations:

1. cos ( ϕ + 1 ) = sin ( ϕ + p ) cos ( ϕ + 1 ) = sin ( ϕ + p )
2. sin ( ϕ ) = - sin ( ϕ + p ) sin ( ϕ ) = - sin ( ϕ + p )
3. tan ( ϕ ) = sin ( ϕ ) sin ϕ + p tan ( ϕ ) = sin ( ϕ ) sin ϕ + p
4. cos ( ϕ ) = - tan ( ϕ + p ) · cos ( ϕ - p ) cos ( ϕ ) = - tan ( ϕ + p ) · cos ( ϕ - p )

(tan(ϕ)tan(ϕ) is defined in exercise 1)

## Trapezoid rule for estimating area (Integration)

### Basic Method

A useful tool for analyzing curves is finding the area underneath them. When we have an unknown combination of waves, we can estimate the area under the curve using the trapezoid rule. We will use f(t)=sin(t)f(t)=sin(t) as an example. If we were to estimate part of the area under the curve with one trapezoid, we might do the following:

We have labeled the two heights, h1h1 and h2h2, and the length of the base bb. The area of the square is:

area of square = ( base ) × ( height ) = b · h 1 area of square = ( base ) × ( height ) = b · h 1
(7)

The area of the top triangle is:

area of triangle = ( base ) × ( height ) 2 = b · ( h 2 - h 1 ) 2 area of triangle = ( base ) × ( height ) 2 = b · ( h 2 - h 1 ) 2
(8)

The total area of the trapezoid is then:

area of trapezoid = b · h 1 + b · ( h 2 - h 1 ) 2 = b · ( h 2 + h 1 ) 2 area of trapezoid = b · h 1 + b · ( h 2 - h 1 ) 2 = b · ( h 2 + h 1 ) 2
(9)

If we know that the two points on the xx-axis are t1t1 and t2t2, then b=t2-t1b=t2-t1. In the figure above, t1=.5t1=.5 and t2=1.5t2=1.5. Then the heights follow from the function: h1=f(t1)=sin(t1)h1=f(t1)=sin(t1) and h2=f(t2)=sin(t2)h2=f(t2)=sin(t2). Thus in general, the area of a trapezoid approximating the area under ff between the points t1t1 and t2t2 is:

area of trapezoid = b · ( h 2 + h 1 ) 2 = ( t 2 - t 1 ) f ( t 1 ) + f ( t 2 ) 2 area of trapezoid = b · ( h 2 + h 1 ) 2 = ( t 2 - t 1 ) f ( t 1 ) + f ( t 2 ) 2
(10)

In order to get a good estimate, we split up the domain of the function f(t)f(t) into several intervals [ti,ti+1][ti,ti+1]. For each interval, we calculate the area of the trapezoid that approximates the area under that curve. For example, we could approximate f(t)f(t) over [0,1][0,1] using four equal intervals. This would look like:

In this case, our estimate would be

approximate area of curve = ( t 2 - t 1 ) f ( t 1 ) + f ( t 2 ) 2 + + ( t 4 - t 3 ) f ( t 3 ) + f ( t 4 ) 2 approximate area of curve = ( t 2 - t 1 ) f ( t 1 ) + f ( t 2 ) 2 + + ( t 4 - t 3 ) f ( t 3 ) + f ( t 4 ) 2
(11)

As we take smaller and smaller intervals, our approximation will get better, because there will be less space between the trapezoids and the curve. We can prove that the trapezoid rule given order 2' convergence–that is, if we cut our intervals in half, our error gets four times smaller.

In the general case, if we split up the domain of the function at points {t1,t2,,tn}{t1,t2,,tn}, then the rule for the estimate is

approximate area of curve = ( t 1 - t 2 ) f ( t 1 ) + f ( t 2 ) 2 + + ( t n - t n - 1 ) f ( t n - 1 ) + f ( t n ) 2 approximate area of curve = ( t 1 - t 2 ) f ( t 1 ) + f ( t 2 ) 2 + + ( t n - t n - 1 ) f ( t n - 1 ) + f ( t n ) 2
(12)

This formula can be further reduced, which is the subject of Exercise 2.1.

### Coding the trapezoid rule

Here we present a code that uses the trapezoid rule to find the area under any function we provide. We need a vector x that holds the values of the domain, for example x = 0:.01:pi. We then need a vector y that holds the function values at those x points, for example y = sin(x).

function curve_area = mytrapz(x, y, fast)
% function curve_area = mytrapz(x, y, fast)
%
% mytrapz.m performs the trapezoid rule on the vector given by x and y.
% Input:
%   x - a vector containing the domain of the function
%   y - a vector containing values of the function corresponding to the

curve_area = 0;

%loop through and add up trapezoids for as many points as we are given
for n = 2 : numel(x)


We start the code with zero area under the curve, since we haven't counted anything yet. Then we create a for loop to count each triangle individually. As we see above, more trapezoids leads to better answers, so we want to use as many trapezoids as we possibly can. In this situation, that means using every point in x and y. The function numel simply counts the number of elements in x. We then calculate the area of the current triangle (within the loop):

    height = (y(n) + y(n-1))/2;     %average height of function across interval
base = x(n) - x(n-1);      	    %length of interval


The area, as in Equation Equation 10, is base times height. The height is the average of the height of the function at the two corner points. The length of the interval is the difference between the two corner points.

Lastly, we compute the area of the trapezoid and add it to our answer:

   trap_area = base * height; 	    %area of trapezoid
curve_area = curve_area + trap_area;    %add to continuing sum
end


The end statement closes the loop over triangles.

To check the accuracy of our code, we test it with many examples, one of which is shown here. The function trapz() is a built-in function that performs the same operation:

>> x = 0:.01:pi;
>> y = sin(x);
>> mytrapz(x,y)

ans =  1.99998206504367

>> trapz(x,y)

ans =  1.99998206504367


Additionally, we know that 0πsin(x)dx=20πsin(x)dx=2, so the answer we are getting is very close to the true value.

A simpler implementation of the code is included at the end, if we specify a third argument fast':

    %alternate (fast) implementation
curve_area = sum(y(2:end-1).*(x(3:end) - x(1:end-2)));
curve_area = curve_area + y(1)*(x(2) - x(1)) + y(end)*(x(end) - x(end-1));
curve_area = curve_area/2;


The explanation for this code is left as exercise 2.1.

### Integration and the Trapezoid Rule

The trapezoid rule is a way of estimating the area under a function. This is exactly what we call an integral. The integral of f(x)f(x) from aa to bb looks like this:

a b f ( x ) d x a b f ( x ) d x
(13)

For well-behaved functions, evaluation of the trapezoid rule approaches the true value of the integral (true area underneath the curve) as we evaluate smaller and smaller intervals. Of course, there are other ways to characterize integration, namely the Fundamental Theorem of Calculus. For more on this, see the Connexions module Integration, Average Behavior: The Fundamental Theorem of Calculus.

### Exercises

2.1 Show that Equation Equation 12 is equivalent to the “fast" version of the code:

approximate area of curve = i = 2 n - 1 f ( t 1 ) ( x i + 1 - x i - 1 ) + f ( t 1 ) ( x 2 - x 1 ) + f ( t n ) ( x n - x n - 1 ) approximate area of curve = i = 2 n - 1 f ( t 1 ) ( x i + 1 - x i - 1 ) + f ( t 1 ) ( x 2 - x 1 ) + f ( t n ) ( x n - x n - 1 )
(14)

2.2 Test the accuracy of mytrapz.m for y=sin(x)y=sin(x) (or some other function) as you increase the number of trapezoids. Plot your result on a loglog() plot: on the xx-axis, number of trapezoids, and on the yy-axis, error.

## Fourier Method for decomposing signals

We now introduce the key mathematical ideas that will allow us to break down signals into their component frequencies.

### Building the Tools

Suppose the function f(t)f(t) consists of several sine waves added together:

f ( t ) = a 1 sin ( 1 · 2 π t ) + a 2 sin ( 2 · 2 π t ) + + a n sin ( n · 2 π t ) f ( t ) = a 1 sin ( 1 · 2 π t ) + a 2 sin ( 2 · 2 π t ) + + a n sin ( n · 2 π t )
(15)

Suppose we know f(t)f(t) over some interval, say [0,1][0,1], and we want to find ajaj. Before we jump into this, we need to calculate the value of a certain integral. For integers hh and kk, with hkhk,

0 1 sin ( h · 2 π t ) sin ( k · 2 π t ) d t = 0 1 1 2 cos ( ( h - k ) · 2 π t ) - cos ( ( h + k ) · 2 π t ) d t trigononmic identity (substitution) = 1 2 ( h - k ) 2 π t sin ( ( h - k ) · 2 π t ) - 1 2 ( h + k ) 2 π sin ( ( h + k ) · 2 π t ) | t = 0 1 integral of cosine is sine = 0 evaluate from 0 to 1 0 1 sin ( h · 2 π t ) sin ( k · 2 π t ) d t = 0 1 1 2 cos ( ( h - k ) · 2 π t ) - cos ( ( h + k ) · 2 π t ) d t trigononmic identity (substitution) = 1 2 ( h - k ) 2 π t sin ( ( h - k ) · 2 π t ) - 1 2 ( h + k ) 2 π sin ( ( h + k ) · 2 π t ) | t = 0 1 integral of cosine is sine = 0 evaluate from 0 to 1
(16)

0 1 sin ( h · 2 π t ) 2 d t = 0 1 1 2 ( 1 - cos ( 2 h · 2 π t ) ) d t trigononmic identity (substitution) = x 2 - 1 ( 2 h · 2 π ) sin ( 2 h · 2 π t ) | t = 0 1 integral of cosine is sine = 1 / 2 evaluate from 0 to 1 0 1 sin ( h · 2 π t ) 2 d t = 0 1 1 2 ( 1 - cos ( 2 h · 2 π t ) ) d t trigononmic identity (substitution) = x 2 - 1 ( 2 h · 2 π ) sin ( 2 h · 2 π t ) | t = 0 1 integral of cosine is sine = 1 / 2 evaluate from 0 to 1
(17)

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

0 1 f ( t ) sin ( j · 2 π t ) d t = 0 1 a 1 sin ( 1 · 2 π t ) sin ( j · 2 π t ) + + a j sin ( j · 2 π t ) sin ( j · 2 π t ) + + a n sin ( n · 2 π t ) sin ( j · 2 π t ) d t ( substitute for f and distribute ) = 0 1 a 1 sin ( 1 · 2 π t ) sin ( j · 2 π t ) d t + + 0 1 a j sin ( j · 2 π t ) sin ( j · 2 π t ) d t + + 0 1 a n sin ( n · 2 π t ) sin ( j · 2 π t ) d t ( the integral of a sum is the sum of integrals ) = 0 + + a j / 2 + + 0 ( all integrals are zero except for the j th integral by the above equations ) = a j / 2 0 1 f ( t ) sin ( j · 2 π t ) d t = 0 1 a 1 sin ( 1 · 2 π t ) sin ( j · 2 π t ) + + a j sin ( j · 2 π t ) sin ( j · 2 π t ) + + a n sin ( n · 2 π t ) sin ( j · 2 π t ) d t ( substitute for f and distribute ) = 0 1 a 1 sin ( 1 · 2 π t ) sin ( j · 2 π t ) d t + + 0 1 a j sin ( j · 2 π t ) sin ( j · 2 π t ) d t + + 0 1 a n sin ( n · 2 π t ) sin ( j · 2 π t ) d t ( the integral of a sum is the sum of integrals ) = 0 + + a j / 2 + + 0 ( all integrals are zero except for the j th integral by the above equations ) = a j / 2
(18)

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

a j = 2 0 1 f ( t ) sin ( j · 2 π t ) d t a j = 2 0 1 f ( t ) sin ( j · 2 π t ) d t
(19)

This derivation is valid for as large of nn 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=b1cos(1·2πt)+b2cos(2·2πt)+f=b1cos(1·2πt)+b2cos(2·2πt)+, then we can recover the coefficients in the same way (Exercise 3.1):

b j = 2 0 1 f ( t ) cos ( j · 2 π t ) d t b j = 2 0 1 f ( t ) cos ( j · 2 π t ) d t
(20)

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

f = sin ( 2 π t ) + sin ( 2 · 2 π t ) + + cos ( 2 π t ) + cos ( 2 · 2 π t ) + f = sin ( 2 π t ) + sin ( 2 · 2 π t ) + + cos ( 2 π t ) + cos ( 2 · 2 π t ) +
(21)

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 ff 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:

T = 5;% duration of signal
dt = 0.001;     % time between signal samples
t = 0:dt:T;
N = length(t);
y = 2.5*sin(3*2*pi*t) - 4.2*sin(4*2*pi*t);    % a 2-piece wave
plot(t,y)
xlabel('time   (seconds)')
ylabel('signal')

for f = 1:5,     % compute the amplitudes as ratios of areas
a(f) = trapz(t,y.*sin(f*2*pi*t))/trapz(t,sin(f*2*pi*t).^2);
end


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(t)cos(j2πt)f(t)cos(j2πt), as in equation Equation 19. 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.

figure
plot(1:5,a,'ko')    % plot the amplitudes vs frequency
hold on
plot(1:5, [0 0 2.5 -4.2 0], 'b*')


As we can see from Figure 7, our method works well.

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.

### Estimating other functions

You may be pleasantly surprised to hear that any periodic function can be broken down into an infinite sum of sine and cosine waves. The elements of such a sum are together called a Fourier Series. The method for recovering the coefficients of the Fourier Series for any periodic function are the same as in equation Equation 19. Thus we create a function that will find the decomposition of any function. We do this in myfourier.m. The entire code is in the appendix, but it uses very similar ideas as in myfreq() (above).

function [freq mag] = myfourier(y, sr, use_fft, plot_on)
(...some code skipped...)

for n = 1 : numel(freq)
%obtain coefficient for freqency nth frequency [freq(n)]
sinmag(n) = mytrapz(t, y.*sin(freq(n)*2*pi*t), 1);
cosmag(n) = mytrapz(t, y.*cos(freq(n)*2*pi*t), 1);
end


As a test, we compare the “answer" fft code with our own for an arbitrary function and find good agreement:

### Exercises

3.1 Suppose f=b1cos(1*2πt)+b2cos(2*2πt)+f=b1cos(1*2πt)+b2cos(2*2πt)+. Then prove Equation Equation 20 that helps recover the coefficients:

b j = 2 0 1 f ( t ) cos ( j · 2 π t ) d t b j = 2 0 1 f ( t ) cos ( j · 2 π t ) d t
(22)

3.2 Suppose ff is an infinite sum of sine and cosine waves:

f = sin ( 2 π t ) + sin ( 2 · 2 π t ) + + cos ( 2 π t ) + cos ( 2 · 2 π t ) + f = sin ( 2 π t ) + sin ( 2 · 2 π t ) + + cos ( 2 π t ) + cos ( 2 · 2 π t ) +
(23)

Show that Equations Equation 19 and Equation 20, which give expressions for ajaj and bjbj, are still valid.

## Spectrogram Applications to EEG

Here is the application we have been seeking the whole time. Given an EEG signal, we would like to find its corresponding Fourier series. We have just built up the tools to decompose the signal into component frequencies. This is useful, but the brain changes over time, and we want to be able to track changes in the signal over time. One way to do this is the spectrogram, or the short-term Fourier transform. It allows us to watch the way a signal changes over time.

### Short-time Fourier Transform and the Spectrogram

The idea behind the Short-time Fourier Transform is that we break up the signal into several time windows, then take the Fourier transform of each one of these. In Figure 10, we analyze the function f(x)=cos(4·2πt)+3cos(10·2πt)+5t2+5sin(35·2πt)f(x)=cos(4·2πt)+3cos(10·2πt)+5t2+5sin(35·2πt) from time t=0t=0 to t=1t=1. First, we break it up into 10 segments of 100 ms each. Each segment corresponds to one block in the “Time" dimension. The line in that block is the amplitudes of sine and cosine waves from the Fourier transform of the function over that time window. The colors on the plane below it represent the amplitude of that frequency in that time interval. The higher the amplitude, the closer to red, and the lower, the closer to blue. The resulting image that we get is called a spectrogram.

To code something like this, we consider each time window separately. This takes the form of a loop:

function [stft_plot freq tm] = my_stft(y, dt, win_len)
(...some code skipped...)
for n = 1:Nwindow
%isolate the part of the signal we want to deal with
sig_win = y((n-1)*win_len + 1 : n*win_len);


Within this window, we perform the Fourier Transform (using our homemade code!) and record it:

    %perform the fourier transform
[mg freq] = myfourier(sig_win, dt, 1);
sm(:,n) = mg(:,1);
cm(:,n) = mg(:,2);
end


We then plot the results on a 2-D plane using imagesc(). This function associates colors with the values in the matrix, so that you can see where the values in the matrix are larger. Thus we get a pretty picture as above.

### Testing the Spectrogram

For the first test, we will try a simple sin wave.

For the second test, we will use a more interesting function. Below we have plotted sin(6*2πt)+sin(20*2πt)+2*sin(et/1.5)sin(6*2πt)+sin(20*2πt)+2*sin(et/1.5) over [0,10][0,10].

This function is interesting because it contains a frequency component that is changing over time. While we have waves at a constant 6 and 20 Hertz, the third component speeds up as tt gets larger and the exponential curve gets steeper. Thus for the plot we expect to see a frequency component that is increasing. This is exactly what we see in Figure 12–two constant bands of frequency, and one train of frequency that increases with time.

>> dt = 1e-4;
>> t = 0:sr:10;
>> y = sin(6*2*pi*t)+sin(20*2*pi*t)+2*sin(exp(t/1.5));
>> my_stft(y, dt, 5000);


### Application to EEG data

For the final section, we will analyze actual brain waves. We recorded from and EEG, and got the signal in Figure 13.

To analyze, we find the time-step in the data, then call mysgram(). This gives us the plot below.

Compare the spectrogram to the raw signal. What do you notice? Perhaps the most notable change is the significant increase in signal magnitude near 18 seconds. This is reflected in the spectrogram with the brighter colors. Also, several "dominant" frequencies emerge. Two faint bands appear at 10 Hz 4 Hz for the first half of the signal. In the last section, a cluster of frequencies between 6 and 10 Hz dominate. Many of the patterns are hidden behind the subtleties in the data, however. Decoding the spectrogram is at least as difficult and creating it. Indeed, it is much less defined. The next section will explore these rules in the context of an interesting application.

#### Application: Driving a Car

One application of decoding brain waves is giving commands to a machine via brainwaves. To see developing work in this field, see this video of the company NeuroSky. Of the many machines we could command, we choose here to command a virtual car (some assembly required) that goes left or right. As above, the decision rule for such a program is not obvious. As a first try, we can find the strongest frequency for each time section and compare it to a set value. If it is lower, the car moves left, and if higher, the car moves right. The following code implements this rule:


N = numel(eeg_sig);
win_len = 3 * round(N / 60);
n = 0;
freq_criterion = 8;

while (n+3) * win_len / 3  <= N %for each time window
%define the moving window and isolate that piece of signal
sig_win = eeg_sig(round(n * win_len / 3) + (1:win_len));

%perform fourier analysis
[freq raw_amps] = myfourier(sig_win, dt, 1);
%only take positive frequencies
freq = freq((end/2+1):end);
%add sine and cosine entries togethef
amps = abs(sum(raw_amps(end/2:end,1), 2));

%find the maximum one
[a idx] = max(amps);

%find the frequency corresponding to the max amplitude
max_freq = freq(idx);

%decided which way the car should move based on the max frequency
if max_freq < freq_criterion;
fprintf('Moving left.... (fmax = %f)\n', max_freq);
%this is where we put animation code
else
fprintf('Moving right.....(fmax = %f)\n', max_freq);
%this is where we put animation code
end

pause(.5); %for dramatic effect :)
n = n + 1;
end


Another approach would be to find the center of mass of the data. This decision rule is more comprehensive because it takes into account all data from the Fourier transform, not just the maximum value. In order to find the average weight of all amplitudes, we change the inner part of code to the following (starting with "%find the maximum one"):

...
%find the frequency corresponding to the "average" amplitude
avg_freq = sum(freq.*amps)/(sum(amps)*df);

%decided which way the car should move based on the max frequency
if avg_freq < freq_criterion;
...

One can imagine myriad other ways to approach this problem. Many strategies have been developed, but the question is open-ended. A natural next step is for the reader to think of new ways to interpret spectrogram data. The most effective characterizations probably have yet to be discovered!

## Conclusion

In this module we developed the tools to decompose an arbitrary signal, such as an EEG, into is component frequencies. We began with sine waves, established the trapezoid scheme, and finally introduced Fourier analysis. This same flavor of analysis is used in many other settings, too–see the related documents.

## Code

### Code for mytrapz.m

function curve_area = mytrapz(x, y, fast)
% function curve_area = mytrapz(x, y, fast)
%
% mytrapz.m performs the trapezoid rule on the vector given by x and y.
%
% Input:
%   x - a vector containing the domain of the function
%   y - a vector containing values of the function corresponding to the
%   values in 'x'

if nargin < 3
curve_area = 0;

%loop through and add up trapezoids for as many points as we are given
for n = 2 : numel(x)
height = (y(n) + y(n-1))/2;     %average height of function across interval
base = x(n) - x(n-1);           %length of interval
trap_area = base * height;      %area of trapezoid
curve_area = curve_area + trap_area;    %add to continuing sum
end

elseif fast
%alternate (fast) implementation
xvals = x(3:end) - x(1:end-2);
yvals = y(2:end-1);
curve_area = yvals(:)'*xvals(:);
curve_area = curve_area + y(1)*(x(2) - x(1)) + y(end)*(x(end) - x(end-1));
curve_area = curve_area/2;
end


### Code for myfreq.m


% myfreq.m
%
% find the frequencies and amplitudes at which a wave is "vibrating"
%
% Contrast simple (but laborious) trapezoid computations to the fast
% and flexible built-in fft command (fft stands for fast Fourier
% transform).
% To make full sense of this we will need to think about complex
% numbers and the complex exponential function.
%

T = 5;% duration of signal
dt = 0.001;     % time between signal samples
t = 0:dt:T;
N = length(t);
y = 2.5*sin(3*2*pi*t) - 4.2*sin(4*2*pi*t);    % a 2-piece wave
plot(t,y)
xlabel('time   (seconds)')
ylabel('signal')

for f = 1:5,     % compute the amplitudes as ratios of areas
a(f) = trapz(t,y.*sin(f*2*pi*t))/trapz(t,sin(f*2*pi*t).^2);
end

figure
plot(1:5,a,'ko')    % plot the amplitudes vs frequency
hold on
plot(1:5, [0 0 2.5 -4.2 0], 'b*')

figure(34)
f = (0:N-1)/T;% fft frequencies
sc = N*trapz(t,sin(2*pi*t).^2)/T;   % fft scale factor
A = fft(y);
newa = -imag(A)/sc;
plot(f,newa,'r+')

y = y + 3*cos(6*2*pi*t);    % add a cosine piece
figure(1)
hold on
plot(t,y,'g')    % plot it
hold off
legend('2 sines','2 sines and 1 cosine')

figure(2)
A = fft(y);           % take the fft of the new signal
newa = -imag(A)/sc;
plot(f,newa,'gx')
b = real(A)/sc;
plot(f,b,'gx')
xlim([0 7])   % focus in on the low frequencies
hold off
xlabel('frequency  (Hz)')
ylabel('amplitude')
legend('by hand','by fft','with cosine')


### Code for myfourier.m

% function [mag freq] = myfourier(y, dt, use_fft)
%
% myfourier.m decomposes the signal 'y', taken with sample interval dt,
% into its component frequencies.
%
% Input:
%
%   y       -- signal vection
%   dt      -- sample interval (s/sample) of y
%   use_fft -- if designated, use matlab's fft instead of trapezoid method
%
% Output:
%
%   freq -- frequency domain
%   mag  -- magnitude of frequency components of y corresponding to 'freq'

function [freq mag] = myfourier(y, dt, use_fft)

y = y(:);
N = numel(y); %number of samples
T = N*dt;                %total time
t = linspace(0,T,N)';     %reconstruct time vector
half_N = floor(N/2);  %ensures that N/2 is an integer
if mod(N,2) %treat differently if f odd or even
freq = (-half_N:half_N)'/T;   %fft frequencies
else
freq = (-half_N:half_N-1)'/T;   %fft frequencies
end

if nargin < 3  %perform explicit Fourier transform
sinmag = zeros(size(freq));  %vector for component magnitudes
cosmag = zeros(size(freq));  %vector for component magnitudes

%loop through each frequency we will test
for n = 1 : numel(freq)
%obtain coefficient for freqency 'freq(n)'
sinmag(n) = mytrapz(t, y.*sin(freq(n)*2*pi*t), 1);
cosmag(n) = mytrapz(t, y.*cos(freq(n)*2*pi*t), 1);
end

%scale to account for sample length
scale_factor = mytrapz(t, sin(2*pi*t).^2);
sinmag = sinmag / scale_factor;
cosmag = cosmag / scale_factor;
mag = [sinmag(:)  cosmag(:)];

elseif use_fft   %use built-in MATLAB fft() for speed
fft_scale_factor = mytrapz(t, sin(2*pi*t).^2) * N / T;
A = fft(y);
mag(:,1) = -imag(A)/fft_scale_factor;
mag(:,2) = real(A)/fft_scale_factor;
mag = circshift(mag, half_N);
end


### Code for mysgram.m

%
% function [stft_plot freq tm] = my_stft(y, dt, Nwindow)
%
% my_stft splits the signa 'y' into time windows, the breaks each
% segment into its component frequencies.  See "Short-time Fourier Transform"
%
%
% Input:
%   y       -- signal
%   dt      -- sample interval
%   Nwindow -- number of time intervals to analyze
%
% Output:
%   stft_plot -- values plotted in the spectrogram
%   freq      -- frequency domain
%   tm        -- time domain

function [stft_plot freq tm hh] = mysgram(y, dt, Nwindow)

%count the number of windows
N = numel(y);
win_len = floor(N/Nwindow);
sm = zeros(win_len, Nwindow);
cm = zeros(win_len, Nwindow);
tm = linspace(0, numel(y) * dt, Nwindow);

%for each window
for n = 1:Nwindow
%isolate the part of the signal we want to deal with
sig_win = y((n-1)*win_len + 1 : n*win_len);
%perform the fourier transform
[freq mg] = myfourier(sig_win, dt, 1);
sm(:,n) = mg(1:win_len,1);
cm(:,n) = mg(1:win_len,2);
end

stft_plot = abs(sm + cm);
stft_plot = stft_plot(end/2:end, :);

%plot the fourier transform over time
hh = imagesc(tm, freq(round(end/2):end), stft_plot);
title('Spectrogram', 'FontSize', 20)
xlabel('time', 'FontSize', 16)
ylabel('frequency', 'FontSize', 16)
set(gca, 'ydir', 'normal')

%just look at lower frequencies
ylim([0-win_len/2 50+win_len/2])


## Content actions

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

#### Collection to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks

#### Module to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks