% ### EXautocorr4.m ###     12.02.14
% Demonstrates how to use autocorrelation to ascertain 'how long does it
% take for a noisy signal to repeat itself?'
% **NOTE** Code not yet complete (need some way to extract the spacing of
% the appropriate maxima from the autocorrelation sans a FFT)
clear
% ---
maxN= 8192;     % largest possible buffer length
SR= 44100;      % sample rate [Hz]
scaleN= 2.05;    % scale factor for noise
R= 5;          % number of repeats
% -------------------------------------------------
N= ceil(rand(1)*maxN);
x= linspace(0,N,N); % determine range
t= x/SR;        % convert to a time
polyS= ceil(10*rand(1));   % randomly determine polynomial order
Pv= 2*rand(polyS,1)-1;  % polynomial coefficients
freqR= 100*(1+scaleN*randn(1));     % randomize up a bit the sin freq. (but not from run to run)
dt= 1/SR;  % spacing of time steps
nPts= N*R;      % total # of points
freq= SR*(0:nPts/2)./nPts;    % create a freq. array (for FFT bin labeling)
% ---
% kludge way to build up signal; tamps down at edges (then adds noise)
y= []; time= [];
for nn=1:R
    temp= hanning(N)'.*(polyval(Pv,t)+sin(freqR*2*pi*t))+scaleN*randn(N,1)';
    y= [y temp]; % pseudo-random 'data'
    
end
% ---
% autocorrelation (via Matlab's cross-correlation routine or a home-built function; 
% the two yhield the same answer, but Matlab's function is faster)
if (1==1),  hh= xcorr(y,y);    
else hh= correlate1(y,y); end
thh= linspace(-N*R/SR,N*R/SR,numel(hh));
hh= hh/max(hh(:));  % normalize to unity
% +++
% plot time waveform
figure(1); clf; plot(linspace(0,numel(y)/SR,numel(y)),y);  hold on; grid on;
xlabel('Time [s]'); ylabel('Signal'); title('Original waveform');
% +++
% plot autocorrelation
figure(2); clf; plot(thh,hh); hold on; grid on; xlabel('Lag [s]'); ylabel('AC'); title('Autocorrelation');
% ---
% now we need some means to extracts the spacing between the relevant maxima...
tf= numel(hh);  % kludge to handle indexing 'end'?
aCorr= hh(ceil(tf/2):tf);    % toss half (redundant due to symmetry)
specAC= rfft(aCorr); specW= rfft(y);
% +++
% plot spectra for comp.
figure(3); clf;
subplot(211); plot(freq,db(abs(specW).^2)); hold on; grid on; title('Power spectrum')
xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]');
subplot(212); plot(freq,db(specAC)); hold on; grid on; title('FT of autocorrelation')
xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]');