format short
format compact
echo on
% A.1- Having saved audio1.mat in an accessible folder, we load it:
% MATLAB version:
% load('audio1');
% octave version:
load -force 'audio1.mat'
% A.2- Now we plot the vector S:
plot(S,'b')
pause
% A.3- Given that S is sampled at 8 KHz, 
% with the first sample representing time 1/8000 
% and the 8000th sample representing time 1, 
% create a second vector of the samples 
% spanning the time period from .275 to .325 seconds inclusive.
% n/8000 = .275, therefore n = .275*8000, etc.
S1 = S(.275*8000:.325*8000);
% Plot it. 
plot(S1,'b');
pause
% A.4- Create a third vector corresponding to the 10399 samples of S
% scaled so that the minimum value is -1 and the maximum value is 1.
S2 = 2*(S-min(S))/(max(S)-min(S)) - 1;
plot(S2,'b');
pause
% or preserving zero, mapping max to 1, and letting min fall where it may:
S3 = S/max(S);
plot(S3,'b');
pause
% A.5- Calculate the serial cross-correlation of 80-sample pieces of S,
% starting at sample 2200, for lags between 1 and 150 samples inclusive.
n=2200;
A=zeros(1,150);
% Note: we use 'ii' rather than 'i' for the index
% because 'i' sometimes gets treated as sqrt(-1)
% and Octave would do so in the expression 'A(i)' below...
for ii=1:150,
    A(ii) = (S(n:(n+79))'*S((n+ii):(n+ii+79))) / ...
        sqrt( (S(n:(n+79))'*S(n:(n+79))) * S(n:(n+ii+79))'*S(n:(n+ii+79)) ) ;
end
plot(A,'b');
pause
% A.6- What is the maximum value of vector A? For what index value does it occur? 
[val,ind] = max(A(2:150))
% now let's show the lag interval corresponding to the maximum
% superimposed on the original waveform of the signal region we worked on:
plot(1:151,S(2200:2350),'b');
hold on
plot([1 151], [0 0], 'k');
[Sval,Sind] = max(S(2200:2230));
plot([Sind Sind], [-20000 25000],'g');
plot([Sind+ind Sind+ind], [-20000 25000],'g');
hold off;
pause
% Just for fun, let's play the sound clip:
% (doesn't work in Octave...)
% soundsc(S);
pause
% A.7- Plot A three times, first with the X axis in samples;
% second with the X axis in (time lag in) seconds (i.e. samples/8000);
% and third with the X axis in (time lag in) Hertz (the inverse of the measure in seconds).
% In the last case, limit the range of values shown to 500 Hz. and below.
  % (note -- Windows version of Octave does not support 'figure()' command)
% figure(1)
plot(A)
title('Serial Cross-correlation: X in samples');
pause
% figure(2)
plot((1:150)/8000,A)
title('Serial Cross-correlation: X in seconds');
pause
% figure(3)
xmax = round(8000/500)
plot(8000./(xmax:150), A(xmax:150))
title('Serial Cross-correlation: X in Hz.');
pause
% Now show the X axis in Hz on a log scale:
% figure(4)
semilogx(8000./(xmax:150),A(xmax:150))
% axis("tight") doesn't seem to work in octave for this plot
axis("tight")
pause
% [Extra Credit] 
% Calculate a version of A based only on the numerator of the equation given 
% (i.e. based only on the inner products of an 80-element stretch of S
% with 150 other 80-element stretches of S.)
A1=zeros(1,150);
for ii=1:150,
    A1(ii) = (S(n:(n+79))'*S((n+ii):(n+ii+79)));
end
% plotyy also doesn't exist in octave...
% figure(1)
% plotyy(1:150,A,1:150,A1);
% a poor substitute -- what would be better?
plot(1:150,(10^10) * A,'g')
title("")
hold on
plot(1:150,A1,'k')
hold off
pause
% Is the maximum at a different place? 
% The plot suggests not --
[val1, ind1] = max(A(2:150))
[val2, ind2] = max(A1(2:150))
%
% NO -- there is a scaling difference of 10^10 or so,
%  but the maxima are at the same lag
%  (as will usually be true since the denominator of the equation
%   is not very strongly dependent on the lag...
A2=zeros(1,150);
for ii=1:150,
        A2(ii)= sqrt( (S(n:(n+79))'*S(n:(n+79))) * S(n:(n+ii+79))'*S(n:(n+ii+79)) ) ;
end
% figure(2)
% plotyy(1:150,A1,1:150,A2);
plot(1:150,A1,'g')
hold on
plot(1:150,A2,'k')
hold off