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