program cross
* Program to generate cross correlation with and with out
* noise
integer*4 mxn
parameter ( mxn = 4096) ! Maximum number of values in TS
real*4 x(mxn), y(mxn), nx(mxn), ny(mxn), rxx(mxn), rxy(mxn)
real*4 ss, sn, c, snr, xs
integer*4 i,j,k, n, ns,ne, num
character*1 tb ! tab character of tab-delimited file.
logical clipping
* Set the sample size (n) and where we will start correlation (ns)
n = 1024
ns = 256
C n = 4096
C ns = 256
* Set sample number of end of correlation (symmetric) and compute
* the number of samples
ne = n-ns
num = ne - ns
* Set clipping true to make samples be +-1, set false to keep original
* samples
clipping = .true.
* Set the tab character
tb = char(9)
* Set the SNR desired and compute the sigmas of the signal and
* noise to acheive this SNR (note variance of uniform distribution
* between -0.5 and 0.5 is 1/12.
SNR = 0.1
if( SNR.lt.1000 ) then
c = 12.d0/(1+SNR)
else ! Take to be infinite SNR
c = 0.d0
endif
* Example values
c c = 10.90909091d0 ! SNR = 0.1
c c = 6.d0 ! SNR = 1
* Compute standard deviation of signal (ss) and of noise (sn)
* (Values choosen so that sum of signals will have variance of 1.0)
ss = sqrt(12.d0-c)
sn = sqrt(c)
* Generate x and y time series. Y is simple x lagged by 50 steps.
* Each time series has independent noise.
* Data is written to unit 20, the auto and cross correlations are
* written to unit 21.
if( clipping ) then
write(20,'(5a)') 'Sample',tb,'X Clipped',tb,'Y Clipped'
else
write(20,'(5a)') 'Sample',tb,'X',tb,'Y'
endif
* Loop over n-samples generating the time series
do i = 1, n
* Value of signal at this time
xs = (rand()-0.5d0)*ss
* Add the noise to the signal
x(i) = xs + (rand()-0.5d0)*sn
* generate y as a lagged version of x with independent noise
if( i+50.le.n ) then
y(i+50) = xs + (rand()-0.5d0)*sn
endif
* Generate random values for y before the start of the lagged
* portion (x also has 50 samples at the end that do not appear
* in the y-timeseries.
if( i.lt.50) then
y(i) = (rand()-0.5d0)*ss + (rand()-0.5d0)*sn
endif
* If we are clipping replace values with +-1 depending on sign
if( clipping ) then
if ( x(i).ge.0 ) then
x(i) = 1.0
else
x(i) = -1.0
end if
if ( y(i).ge.0 ) then
y(i) = 1.0
else
y(i) = -1.0
end if
endif
* Write out the data for plotting (offset values by +-2 so that
* separate time series can be seen).
write(20,100) i, tb, x(i)+2.0,tb, y(i)-2.0
100 format(i5,a,F10.6,a,f10.6)
end do
**** Now form the correlation function
* Write out the header line for the output to unit 21
if( sn.gt.0 ) then
if( clipping ) then
write(21,'(5a,f6.2)')
. 'Lag',tb,'\\fr\\nxx Clip',tb,
. '\\fr\\nxy Clip SNR ',snr
else
write(21,'(5a,f6.2)')
. 'Lag',tb,'\\fr\\nxx',tb,'\\fr\\nxy SNR ',snr
endif
else
if( clipping ) then
write(21,'(5a)')
. 'Lag',tb,'\\fr\\nxx Clip',tb,
. '\\fr\\nxy Clip SNR Infinite'
else
write(21,'(5a)')
. 'Lag',tb,'\\fr\\nxx',tb,'\\fr\\nxy SNR Infinite'
endif
end if
* Now do the correlation just over +-100 laggs about the center
* (Offset is 50 so will appear half way to one side).
do k = -100,100
rxx(k+101) = 0.d0
rxy(k+101) = 0.d0
do j = ns,ne
* Form autocorrelation of x and cross correlation of x and y
rxx(k+101) = rxx(k+101)+x(j)*x(j+k)
rxy(k+101) = rxy(k+101)+x(j)*y(j+k)
end do
* Nomralize by number of samples
rxx(k+101) = rxx(k+101)/num
rxy(k+101) = rxy(k+101)/num
* Write out the values`
write(21,210) k,tb,rxx(k+101),tb,rxy(k+101)
210 format(i5,a,F10.6,a,f10.6)
end do
end