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=[];
%start doing what calCSD does
%This normalization makes fc as if from (0,1) data.
nw=2/norm(wind)^2;
%Zero padding is done by windowAndFFT
zpad=zeros(numPoints,1);
%Prepare for coarseGrain.
numFreqs=1/2*resampleRate/deltaF;
cz=zeros(numFreqs,1);
jsum=2*deltaF*segmentDuration;
for in=1:numsegments
n3a=randn(numPoints/2,1);
n1a=randn(numPoints/2,1);
n2a=randn(numPoints/2,1);
n4a=randn(numPoints/2,1);
n3b=randn(numPoints/2,1);
n1b=randn(numPoints/2,1);
n2b=randn(numPoints/2,1);
n4b=randn(numPoints/2,1);
na=[wind.*[n1a;n2a];zpad];
nb=[wind.*[n1b;n2b];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
c=c+fc((jn-1)*numFreqs+1:jn*numFreqs);
end
d=[d;c];
n1a=n1b;
n2a=n2b;
[Ta31(:,in),Fa31(:,in)]=psd([n3a;n1a],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag); 
[Ta24(:,in),Fa24(:,in)]=psd([n2a;n4a],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
[Ta12(:,in),Fa12(:,in)]=psd([n1a;n2a],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
[Tb31(:,in),Fb31(:,in)]=psd([n3b;n1b],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag); 
[Tb24(:,in),Fb24(:,in)]=psd([n2b;n4b],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
[Tb12(:,in),Fb12(:,in)]=psd([n1b;n2b],psdFFTLength,1/deltaT,psdWindow,psdOverlapLength,detrendFlag);
end
Tal=[(Ta31(1,:)+Ta31(end,:))/2 reshape(Ta31(2:end-1,:),1,(size(Ta31,1)-2)*size(Ta31,2))];
Tar=[(Ta24(1,:)+Ta24(end,:))/2 reshape(Ta24(2:end-1,:),1,(size(Ta24,1)-2)*size(Ta24,2))];
Tac=[(Ta12(1,:)+Ta12(end,:))/2 reshape(Ta12(2:end-1,:),1,(size(Ta12,1)-2)*size(Ta12,2))];
Ta=(Tal+Tar)/2;
Tbl=[(Tb31(1,:)+Tb31(end,:))/2 reshape(Tb31(2:end-1,:),1,(size(Tb31,1)-2)*size(Tb31,2))];
Tbr=[(Tb24(1,:)+Tb24(end,:))/2 reshape(Tb24(2:end-1,:),1,(size(Tb24,1)-2)*size(Tb24,2))];
Tbc=[(Tb12(1,:)+Tb12(end,:))/2 reshape(Tb12(2:end-1,:),1,(size(Tb12,1)-2)*size(Tb12,2))];
Tb=(Tbl+Tbr)/2;
T34=sqrt(Ta'.*Tb');
T12=sqrt(Tac.*Tbc);
%T24=sqrt(Tar.*Tbr);
%Note: normalization here differs from coarseGrain output.
dd=[real(d);imag(d)];
%snr=real(d)./T34);
save('snrdata.mat','d','T34');