Table of Contents
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 ϕ of the triangle. As you can see, if ϕ is the angle of a right triangle with hypotenuse 1 (illustrated by the circle) , sin ( ϕ ) cos ( ϕ )
sincirc.m (copied below) to see how the sine and cosine values relate to the angle
|
% 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:
Aside: The wave described by sin ( t ) 1 / 2 π sin ( 2 π t )
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:
A wave with wavelength ℓ is written as sin ( 2 π t ℓ )
2.
Amplitude measures how high the wave is. We can make a wave of amplitude a by writing:
Sometimes a will be negative. In this case, we still say that the wave has amplitude a , 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 p , we write:
Since a sin wave repeats every 2 π ϕ ):
Aside: A cosine wave is a sine wave shifted back by a π / 2
Every sine (cosine) wave can be described completely by these characteristics. These are shown in Figure 2 (phase = 0
 |
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.
Adding Sine Waves
We can add together multiple sin waves to accomplish different shapes. For example, if we add the wave [ 4 · sin ( 2 π t ) ] [ sin ( 6 · ( 2 π t ) ) ]
 |
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 1 6
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 ) sin ( ϕ ) = sin ( ϕ + 2 π ) p in the following equations by making the above substitutions. Check your answer visually against the figures of sine and cosine waves:
Equations:
(tan ( ϕ )
|
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 )
|
We have labeled the two heights, h 1 h 2 b . 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 x -axis are t 1 t 2 b = t 2 - t 1 t 1 = . 5 t 2 = 1 . 5 h 1 = f ( t 1 ) = sin ( t 1 ) h 2 = f ( t 2 ) = sin ( t 2 ) f between the points t 1 t 2
In order to get a good estimate, we split up the domain of the function f ( t ) [ t i , t i + 1 ] f ( t ) [ 0 , 1 ]
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 { t 1 , t 2 , ⋯ , t n }
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 ) d x = 2
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 ) a to b 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.
Exercises
2.1 Show that Equation Equation 12 is equivalent to the “fast" version of the code:
2.2 Test the accuracy of y = sin ( x ) x -axis, number of trapezoids, and on the y -axis, error.
mytrapz.m for loglog() plot: on the 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 )
Suppose we know f ( t ) [ 0 , 1 ] a j h and k , withh ≠ k
If instead h = k
Now we wish to find the coefficients of f ( t ) f ( t ) sin ( j · 2 π t )
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 π t ) + b 2 cos ( 2 · 2 π t ) + ⋯
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: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 f ( t ) cos ( j 2 π t )
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 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 = b 1 cos ( 1 * 2 π t ) + b 2 cos ( 2 * 2 π t ) + ⋯
3.2 Suppose f is an infinite sum of sine and cosine waves:
Show that Equations Equation 19 and Equation 20, which give expressions for a j b j
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 ) + 3 cos ( 10 · 2 π t ) + 5 t 2 + 5 sin ( 35 · 2 π t ) t = 0 t = 1
 |
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 ( e t / 1 . 5 ) [ 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 t 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:
%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!
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])
Posts
0 comments:
Post a Comment