% COGS501 demo (b)
figure(1); subplot(1,1,1); plot(0);
format short; format compact; clc; echo on;
%
% Now let's calculate a time function of amplitude for the audio file that
% we read just read in.
%
% This will be the "root mean square" amplitude: 
%      Each value will be
%      the square root of the sum of the squared sample values
%      in a short time region of the signal
% We'll use overlapping 20-millisecond analysis windows,
%    advanced 10-milliseconds at a time (or 100 samples per second).
% The window length in samples will be
w = round(.02*SR)
% The length of our signal is
n = length(M)
pause
% How many analysis windows will "fit" in our signal?
% The last sample a window can start on is n-w+1
%     so we want the biggest integer p such that 
%        p*.01*SR <= n-w+1
% It should work to use
p = floor( (n-w+1)/(.01*SR) )
pause
% We can check that this value works
[n round(p*.01*SR)+w]
% and that the next one doesn't:
[n round((p+1)*.01*SR)+w]
pause
% Note that the analysis-window starting points will be (rounded)
% 0 * .01*SR
% 1 * .01*SR
% 2 * .01*SR
% ...
% p * .01*SR
%
% so that there are actually p+1 analyses.
%
pause
home
% OK, how can we do the analysis?
% Well, if we have a vector of w values, e.g.
R = rand(w,1);
% then we can get the sum of squares by doing a simple inner product:
R'*R
%
pause
% Note that this is the same as using an explicit multiply-add loop:
% sumr = 0;
% for(i=1:w)
%    sumr = sumr+R(i)^2;
% end
echo off
sumr = 0;
for(i=1:w)
    sumr = sumr+R(i)^2;
end
echo on
sumr
%
pause
% Of course the mean squared value is just R'*R/w
%    and the square root of that is sqrt(R'*R/w)
%
% And we can use Matlab's "slice" syntax for pulling out
%     the bit of M that we need for each calculation...
pause
% Thus we can calculate a time-function of rms amplitude as follows:
% n = length(M);
% w = round(.02*SR);
% p = floor( (n-w+1)/(.01*SR) );
% aM = zeros(p+1,1);
%
% for(i=0:p)
%    sbase = round(i*.01*SR)+1; srange = sbase:(sbase+w-1);
%    aM(i+1) = sqrt(M(srange)'*M(srange)/w);
% end
%
echo off
n = length(M);
w = round(.02*SR);
p = floor( (n-w+1)/(.01*SR) );
aM = zeros(p+1,1);
%
for(i=0:p)
    sbase = round(i*.01*SR)+1; srange = sbase:(sbase+w-1);
    aM(i+1) = sqrt((M(srange)'*M(srange))/w);
end
echo on
pause
% 
% And plot the result --
% Note that window center times are
% 0*.01 + .5*w/SR
% 1*.01 + .5*w/SR
% 2*.01 + .5*w/SR
% ...
% p*.01 + .5*w/SR
%
plot((0:p)*.01 + .5*w/SR, aM)
%
figure(2); subplot(1,1,1); plot((1:length(M))/SR, M);
pause
% To get amplitude in decibels (DB), we'd use 20*log10(aM):
figure(3); subplot(1,1,1);
plot((0:p)*.01 + .5*w/SR, 20*log10(aM))
pause
%
% Now we can turn what we just did into a function called "getamp"
%     by making a new m-file "getamp.m"
%     containing something like this 
%        (minus the first two chars on each line)
%
%| function y = getamp(X, SR, window, increment)
%| % y = getamp(X, SR, window, increment)
%| %     X -- audio vector
%| %     SR -- sampling rate (in samples/second)
%| %     window -- analysis window (in seconds) (square!)
%| %     increment -- analysis frame increment (in seconds)
%| n = length(X);
%| w = round(window*SR);
%| p = floor( (n-w+1)/(.01*SR) );
%| y = zeros(p+1,1);
%| for(i=0:p)
%|     sbase = round(i*.01*SR)+1; srange = sbase:(sbase+w-1);
%|     y(i+1) = sqrt((X(srange)'*X(srange))/w);
%| end
%
pause
% And we can test it:
aM1 = getamp(M, SR, .02, .01);
figure(4); subplot(1,1,1)
plot((0:p)*.01 + .5*w/SR, aM1)
%
pause
%
home
% Now let's count the number of time in each analysis window
%     that the signal crosses zero: the "zero crossing rate".
% If we multiply each element of a vector by the next element in sequence,
%    the result will be positive if both elements are positive 
%        and also if both elements are negative,
%    but the result will be negative 
%        if the two are on opposite sides of zero.
R = rand(10,1)-.5
R(1:9).*R(2:10)
pause
% A clearer display:
[R [0; R(1:9).*R(2:10)] ]
pause
% If we ask whether the values are less than zero
R(1:9).*R(2:10) < 0
% we get 1's and 0's for True and False,
% so we can count by summing:
sum(R(1:9).*R(2:10) < 0)
%
pause
% So here's a script for calculating zero crossing rates:
%| n = length(M);
%| w = round(.02*SR);
%| p = floor( (n-w+1)/(.01*SR) );
%| zM = zeros(p+1,1);
%| %
%| for(i=0:p)
%|     sbase = round(i*.01*SR)+1; srange = sbase:(sbase+w-1);
%|     zM(i+1) = sum((M(srange).*M(srange+1)) < 0 );
%| end
pause
echo off
n = length(M);
w = round(.02*SR);
p = floor( (n-w+1)/(.01*SR) );
zM = zeros(p+1,1);
%
for(i=0:p)
    sbase = round(i*.01*SR)+1; srange = sbase:(sbase+w-1);
    zM(i+1) = sum((M(srange).*M(srange+1)) < 0 );
end
echo on
% Let's plot it:
figure(5); subplot(1,1,1)
plot((0:p)*.01 + .5*w/SR, zM)
pause
%
% If we look at a shorter stretch, we can see that the amplitude
%    and the zero crossing count often go in opposite directions:
figure(4)
plot((400:600)*.01 + .5*w/SR, aM(400:600))
figure(5)
plot((400:600)*.01 + .5*w/SR, zM(400:600))
pause
% A scatter plot also shows this:
figure(6)
scatter(aM,zM)
pause
home
% From these examples, you've learned how to derive and plot 
%   some simple parameters of a time series,
%   and how to set up a Matlab function for such calculations.