tic
segmentDuration=52;
resampleRate=4096;
deltaF=0.25;
deltaT=1/resampleRate;
numPoints=segmentDuration*resampleRate;
psdFFTLength=resampleRate/deltaF;
psdWindow= hann(psdFFTLength,'periodic');
csdWindow= hann(numPoints,'periodic');
wind=csdWindow;
psdOverlapLength=psdFFTLength/2; 
detrendFlag='none';
%Choose the number of segments to analyze
numsegments=2000;
d=[];
T1L=[];
T2L=[];
T1C=[];
T2C=[];
T1R=[];
T2R=[];
%start doing what calCSD does
n1a=randn(numPoints/2,1);
n1b=randn(numPoints/2,1);
n1c=randn(numPoints/2,1);
n2a=randn(numPoints/2,1);
n2b=randn(numPoints/2,1);
n2c=randn(numPoints/2,1);
[T1l,F]=psd([n1a;n1b],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag); 
[T2l,F]=psd([n2a;n2b],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
[T1c,F]=psd([n1b;n1c],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
[T2c,F]=psd([n2b;n2c],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag); 
%This normalization makes fc as if from (0,1) data, and goes with sqrt(2*jsum) in snr.
%nw=2/norm(wind)^2;
%This normalization accounts for frequency leakage, and goes with 1/sqrt(jsum) in snr.
nw=1/numPoints/deltaF/2;
%Zero padding is done by windowAndFFT
zpad=zeros(numPoints,1);
%Prepare for coarseGrain.
numFreqs=1/2*resampleRate/deltaF;
cz=zeros(numFreqs-1,1);
jsum=2*deltaF*segmentDuration;
w=[1/2;ones(jsum-1,1);1/2];
outputsize=(numFreqs-1)*numsegments;
for in=1:numsegments
na=[wind.*[n1b;n1c];zpad];
nb=[wind.*[n2b;n2c];zpad];
fa=fft(na);
fb=fft(nb);
fa(1)=fa(1)+i*fa(end/2+1);
fb(1)=fb(1)+i*fb(end/2+1);
%This step is from calCSD.
fc=nw*conj(fa(1:end/2)).*fb(1:end/2);
%The normalization reflects the windowing and zero padding.
%Without them, the correct factor would be 4/numPoints.
%Start the coarseGrain calculation.
c=cz;
for jn=1:jsum+1
c=c+w(jn)*fc(jsum/2+jn:jsum:end-3*jsum/2+jn);
end
d=[d;c];
n1d=randn(numPoints/2,1);
n2d=randn(numPoints/2,1);
[T1r,F]=psd([n1c;n1d],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
[T2r,F]=psd([n2c;n2d],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
T1L=[T1L,T1l];
T2L=[T2L,T2l];
T1C=[T1C,T1c];
T2C=[T2C,T2c];
T1R=[T1R,T1r];
T2R=[T2R,T2r];
n1a=n1b;
n2a=n2b;
n1b=n1c;
n2b=n2c;
n1c=n1d;
n2c=n2d;
T1l=T1c;
T1c=T1r;
T2l=T2c;
T2c=T2r;
end
T1L(1,:)=(T1L(1,:)+T1L(end,:))/2;
T1C(1,:)=(T1C(1,:)+T1C(end,:))/2;
T1R(1,:)=(T1R(1,:)+T1R(end,:))/2;
T2L(1,:)=(T2L(1,:)+T2L(end,:))/2;
T2C(1,:)=(T2C(1,:)+T2C(end,:))/2;
T2R(1,:)=(T2R(1,:)+T2R(end,:))/2;
Tal=reshape(T1L(2:end-1,:),outputsize,1);
Tac=reshape(T1C(2:end-1,:),outputsize,1);
Tar=reshape(T1R(2:end-1,:),outputsize,1);
Ta=(Tal+Tar)/2;
Tbl=reshape(T2L(2:end-1,:),outputsize,1);
Tbc=reshape(T2C(2:end-1,:),outputsize,1);
Tbr=reshape(T2R(2:end-1,:),outputsize,1);
Tb=(Tbl+Tbr)/2;
Tadj=sqrt(Ta.*Tb);
Tloc=sqrt(Tac.*Tbc);
%Note: normalization here differs from coarseGrain output.
%dd=[real(d);imag(d)];
snr=real(d)./Tadj/sqrt(jsum);
save('snrdata.mat','d','Tadj');
toc