% create two data sets with random noise with standarad dev 1
% one with two sine waves at period of TTS and TTL, one with just
% a sine wave of period TTS but shifted in time by "delay"
% here the shared sine has a signal to noise of 2:1, while the 
% non-shared sine wave has a 10:1 signal to noise.
% NaN's are placed radomly in the data sets

% a lagged Cross Correlation is performed, and this is Fourier
% transformed to find the Cross Spectrum
% The frequency of the significant peak in the cross
% spectrum (significance determined by Lomb-Scargle)is given 

% set some varaibles for use:
% number of data points
clear;
% set the periods
TTS = 50;  % the shared period
TTL = 120;
% set the number of data points in fake data
npts = 200;
% set the delay of the common sine wave as 10% npts (-1=delay)
delay = -1*npts/10;
% number of NaN values in each data set.
nNaN = 10;
% number of lags for which to calculate the cross correlaiton
nlags = 155;
% set the oversampling factor for the spectral data
ofac=8;
% significance level to search for peaks
siglev = 0.95;

%*****************************************************************
%            create the fake data files
%*****************************************************************

% set time variable
t = (0:npts-1)';

% now the leading variable with two sine waves
leadpos=randn(npts,1);
leadpos = leadpos + 2*sin(2*pi*t/TTS) + 10*sin(2*pi*t/TTL);

%now the data set with single sine wave trailing by 10% of npts 
leadneg=randn(npts,1);
leadneg = leadneg + 2*sin(2*pi*(t+delay)/TTS) ;

% Add nNAN's to the data set at random locations
% give the randum number generator a random seed
rand('twister',sum(100*clock))
% Then set random values of the two data sets to NaN
for jj = 1:nNaN
    index = fix(rand*npts);
    if (index==0) 
        index = 1;
    end
    leadpos(index) = NaN;
        index = fix(rand*npts);
    if (index==0) 
        index = 1;
    end
    leadneg(index) = NaN;
end

% plot the two data sets in figure 1
figure (1);
plot(leadpos); hold, plot(leadneg, 'r'); hold
xlabel('Time Units');

%*****************************************************************
%        Now perform the cross correlation with NaN's
%*****************************************************************

% now perform the cross-correlaiton at nlags and plot the results 
% in figure 2 using the call
% [xcor, lags, upper, lower] = XcorelNaN(leadpos,leadneg,nlags,fig);

% xcor will be the lagged correlation function at 2*nlags + 1
% to account for positive, negative and zero lag, and the routine 
% also returns the + and - 95% confidence levels. If the shaded 
% area of the plot does not include 0, then the value of xcor can
% be taken to be significant at the 95% level.  The value of the 
% correlation can be taken to be xcor +upper -lower at the 95% level

% The location of the first peak in the positive lag direction 
% indicates the number of time units (lags) that the first data in 
% the call, leadpos, leads the second data set in the call, leadneg. 
% In this example, that means that the first maximum in the sine wave 
% common to both data sets reaches a maximum value in leadpos before 
% the sine wave in leadneg.
% The location of the first peak in the negative lag directions shows 
% the number of time units (lags) the second data set leads the first.
% For the values TTS=50 (period of the common sine wave) and delay=-20
% applied to the second data set, we should expect a peak in the cross
% correlation at -30 and +20 lags. That is, I can shift the second 
% sine wave 20 lags ealier in tie to get the first peaks in the sine 
% waves to match up.
% Alternatively, I can shift the second sine wave 30 lags later in time
% in order to have its sine wav's first peak match up with the second
% peak of the first sine wave.  Such is the nature of sines!
fig = 2;
[xcor, lags, upper, lower] = XcorelNaN(leadpos,leadneg,nlags,fig);

%*****************************************************************
%        Now perform the Cross-Spectrum
%*****************************************************************

% Since the cross-correlation is now continuous, we can do an
% FFT that requires no NaN's. This will result in a cross spectrum

% Pad the data set by ofac and go to next highest power of 2.  
% FFT's work with 2^N points, and the zero 
% padding allows higher sampling of the frequency spectrum 
NFFT = 2^nextpow2(ofac*(2*nlags+1));

% set the frequency variable (only plotting 0 to the Nyquist)
f =1/2*linspace(0,1,NFFT/2);

% set a hamming window function to taper the data to 0 and avoid
% side lobes. Remember, the xcor has 2*nlags+1 elements
W=hamming(2*nlags+1);

%perfom the FFT on the windowed, or apodized, data
Y = fft((xcor.*W),NFFT)/(2*nlags+1);

% plot the frequency spectrum in figure 3
% semi-log makes it easier to spot the peak location
figure(3);
semilogx(f,2*abs(Y(1:NFFT/2)))
xlabel('Frequency 1/(Time Units)');
ylabel('Cross Spectral Amplitude');

% for convenience, plot it as a function of period in figure (4)
% semi-log makes it easier to spot the peak location
PER=1./f;
figure(4)
semilogx(PER,2*abs(Y(1:NFFT/2)));
xlabel('Period (Time Units)');
ylabel('Cross Spectral Amplitude');

%*****************************************************************
%        find the freqyency of max in cross-spectrum
%        using Lomb-Scargle to find significant peaks
%*****************************************************************

% to narrow the search range, 
% do a Lomb-Scargle to get significance limits and only look where 
% significance (1-alpha), is above siglev (alpha = 1-siglev) and 
% regions where the amplitude is above 10% of strongest peak
alphatst = 1-siglev;

% call fastlomb((x,t,fig,hifac,ofac)) with windowed correlation function, 
% hifac=1 (0 to nyquist),  ofac set above and figure 5
[LSP,f1,alpha] = fastlomb(xcor.*hamming(2*nlags+1), lags, 5, 1, ofac);

% now get the frrequency range of interest (higher significance than 
% specified and all peaks within 10% of the maximum value of spectrum
f1int=f1((alpha<alphatst) & (LSP>=max(LSP)/10));

% Now go back to the FFT and get the frequencies and 
% FFT complex amplitudes over this FFT frequency range, f
Yint=Y((f >=min(f1int))&(f<=max(f1int)));
fint=f((f >=min(f1int))&(f<=max(f1int)));

% get the value of the frequency of the peak amplitude
freqp=fint(abs(Yint)==max(abs(Yint)));
disp(strcat([' Peak at frequency ',num2str(freqp),' 1/(time-units)']))
% for convenience, convert this to a peiod
perp = 1/freqp;
disp(strcat(['   or ',num2str(perp),' time units',' compared to set value of ',num2str(TTS)]))
% this is the frequency and period of the sine wave common to both data
% sets, even though is has a very low signal to noise.

