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:

Figure 1: Two illustrations of the sine function.

A Sin wave The relation of sine and cosine to a triangle and unit circle (a) (b)

% 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

The sin wave has three primary characteristics:

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

Figure 2: A sin wave.

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);

>>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);

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:

Figure 3: A compound wave

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.

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:

Figure 6: Applying the trapezoid rule to a sine wave.

A single trapezoid. Estimating the area with four trapezoids. (a) (b)

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

The area of the top triangle is:

The total area of the trapezoid is then:

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:

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

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

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

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=2∫0π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.

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:

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.

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

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.

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

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 h≠kh≠k,

If instead h=kh=k,

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:

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

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

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).

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.

Figure 7: Our method gets the same result as Matlab's fft().

Frequencies captured by myfreq(). Frequencies captured by Matlab's fft(). (a) (b)

Figure 8: The function in question.

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.

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:

Figure 9: Comparison of myfft() to true frequencies.

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:

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

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

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.

Figure 10: Visual illustration of 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.

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].

Figure 11: An interesting signal to decompose.

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);

Figure 12: The spectrogram of the above function

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

Figure 13: An EEG wave.

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

Figure 14

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.

  
%load data
load bwave
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!