OpenW # 1
FullW # 1
FontSize = 10
AutoRedraw = 1
PrintScroll = 1
Global Int nd, i, version
Global String b, ok
b = " "
Print "  Numerical Simulation of a Cross-Correlation Test"
Print
Print "  We will compute a signal  s ";
Print " supposed to represent ground motion,"
Print "  and three independent signals nx, ny, nz";
Print " representing instrumental noise."
Print "  Each of them consists of nd";
Print " samples of gaussian white noise."
Input"  ==> Enter the number of samples, nd: ";nd
Global Double sx(nd), sy(nd), sz(nd), nx(nd), ny(nd) , nz(nd), sig(nd)
Global Double gx, gy, gz, rx, ry, rz, pxx, pxy, pxz, pyy, pyz, pzz
Global Double s, z, pxxth, pxyth, pxzth, pyyth, pyzth, pzzth, pss
Global Double nxx, nyy, nzz, w2, twopi, gym, gyp, nxxm, nxxp, nyym, nyyp
'##E#
Global Double psx, psy, psz, nxis, nyis, nzis
version = 2
w2 = Sqrt(2)
twopi = 2. * PI
Print "  We then construct output signals";
Print " sx, sy, sz as sx = gx * s + rx * nx"
Print "  etc. gx is the gain (transfer function)";
Print " of the x channeland nx the rms"
Print "  amplitude of noise added at the output."
Print "  Same for y and z."
Input "  ==> Enter gain factor gx: "; gx
Input "  ==> Enter gain factor gy: "; gy
Input "  ==> Enter gain factor gz: "; gz
Input "  ==> Enter noise amplitude rx: ";rx
Input "  ==> Enter noise amplitude ry: ";ry
Input "  ==> Enter noise amplitude rz: ";rz
Print "  We now apply, in time domain,";
Print " the conventional cross-correlation"
Print "  analysis and the SWT method. For the";
Print " conventional method we"
Print "  would have to estimate the transfer";
Print " function of the y channel a priori."
Print "  The program uses the correct value and two";
Print " values that are plus and"
Print "  minus one percent off."
Print
gyp = gy * 1.01
gym = gy * 0.99
nxis = 0.
nyis = 0.
nzis = 0.
For i = 1 To nd
  s = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd()
  s = s + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() - 6.
  ' s = w2 * Sin(twopi * i / 25)
  z = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd()
  z = z + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() - 6.
  nx(i) = rx * z
  nxis = nxis + nx(i) ^ 2
  z = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd()
  z = z + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() - 6.
  ny(i) = ry * z
  nyis = nyis + ny(i) ^ 2
  z = Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() + Rnd()
  z = z + Rnd() + Rnd() + Rnd() + Rnd() + Rnd() - 6.
  nz(i) = rz * z
  nzis = nzis + nz(i) ^ 2
  sx(i) = gx * s + nx(i)
  sy(i) = gy * s + ny(i)
  sz(i) = gz * s + nz(i)
  sig(i) = s
Next i
pss = 0.
pxxth = gx ^ 2 + rx ^ 2
nxis = nxis / nd
pxx = 0.
pxyth = gx * gy
pxy = 0.
pxzth = gx * gz
pxy = 0.
pyyth = gy ^ 2 + ry ^ 2
nyis = nyis / nd
pyy = 0.
pyzth = gy * gz
pyy = 0.
pzzth = gz ^ 2 + rz ^ 2
nzis = nzis / nd
pzz = 0.
For i = 1 To nd
  pss = pss + sig(i) * sig(i)
  pxx = pxx + sx(i) * sx(i)
  pxy = pxy + sx(i) * sy(i)
  pxz = pxz + sx(i) * sz(i)
  pyy = pyy + sy(i) * sy(i)
  pyz = pyz + sy(i) * sz(i)
  pzz = pzz + sz(i) * sz(i)
Next i
pss = pss / nd
pxx = pxx / nd
pxy = pxy / nd
pxz = pxz / nd
pyy = pyy / nd
pyz = pyz / nd
pzz = pzz / nd
If version = 1 Then
  Print "  White Noise, "; nd; " samples"
  Print "  pxx expected "; pxxth, "actual "; Round(pxx, 3); ""
  Print "  pxy expected "; pxyth, "actual "; Round(pxy, 3); ""
  Print "  pxz expected "; pxzth, "actual "; Round(pxz, 3); ""
  Print "  pyy expected "; pyyth, "actual "; Round(pyy, 3); ""
  Print "  pyz expected "; pyzth, "actual "; Round(pyz, 3); ""
  Print "  pzz expected "; pzzth, "actual "; Round(pzz, 3); ""
Else
  Print "  conventional two-instrument analysis:"
  nxx = pxx - pxy * gx / gy
  nxxm = pxx - pxy * gx / gym
  nxxp = pxx - pxy * gx / gyp
  nyy = pyy - pxy * gy / gx
  nyym = pyy - pxy * gym / gx
  nyyp = pyy - pxy * gyp / gx
  Print "  nxx exact  ", Round(nxis, 3), "nxx estim. "; Round(nxx, 3)
  Print "  gym = gy-1% = ", gym, "nxx estim. "; Round(nxxm, 3)
  Print "  gyp = gy+1% = ", gyp, "nxx estim. "; Round(nxxp, 3)
  Print "  nyy exact  ", Round(nyis, 3), "nyy estim. "; Round(nyy, 3)
  Print "  gym = gy-1% = ", gym, "nyy estim. "; Round(nyym, 3); ""
  Print "  gyp = gy+1% = ", gyp, "nyy estim. "; Round(nyyp, 3); ""
  psx = (pxx - nxx) / gx ^ 2
  psy = (pyy - nyy) / gy ^ 2
  Print "  The power of the coherent signal is "; Round(pss, 3); " and estimated as"
  Print "  from sx: "; Round(psx, 3), "from sy: "; Round(psy, 3), "avg="; Round((psx + psy) / 2, 3)
EndIf
Print
Input"continue?";ok
Print "  Method of Sleeman, Wettum, Trampert:"
nxx = pxx - pxy * pxz / pyz
nyy = pyy - pxy * pyz / pxz
nzz = pzz - pxz * pyz / pxy
Print "  nxx exact  ", Round(nxis, 3), "estim. "; Round(nxx, 3)
Print "  nyy exact  ", Round(nyis, 3), "estim. "; Round(nyy, 3)
Print "  nzz exact  ", Round(nzis, 3), "estim. "; Round(nzz, 3)
Print "   gx exact    ", gx, "estim. "; Round(Sqrt(pxx - nxx), 3)
Print "   gy exact    ", gy, "estim. "; Round(Sqrt(pyy - nyy), 3)
Print "   gz exact    ", gz, "estim. "; Round(Sqrt(pzz - nzz), 3)
psx = (pxx - nxx) / gx ^ 2
psy = (pyy - nyy) / gy ^ 2
psz = (pzz - nzz) / gz ^ 2
Print "  The power of the coherent signal is "; Round(pss, 3); " and estimated as"
Print "  from sx: "; Round(psx, 3), " from sy: "; Round(psy, 3), " from sz: "; Round(psz, 3), "avg="; Round((psx + psy + psz) / 3, 3)
Print
Input "finished?";ok
$DatFile
$ExeName = C:\Programme\GB32\sleeman_ia.Exe
$FileDescription = Beschreibung
$Comments = Autor: Nicht registriert
$CompanyName = Nicht registriert
$LegalCopyright = © 2005 Nicht registriert
$LanguageId = 0407
$FileVersion = %d.%d.%d.%d
$ProductVersion = %d.%d.%d.%d
$DatFile