% SIMULATEH1 Simulated noise at LIGO design sensitivity
%
% SIMULATEH1 produces stationary Gaussian noise with an amplitude spectral
% density approximating the LIGO design sensitivity.  Stationary white
% noise is filtered by a number of shaping filters to approximate the shot
% noise, thermal noise, and seimic noise of the detector.
%
% SIMULATEH1 produces a continuous seamless stream of data which is then
% written to frame files.
%
% The first frame file contains startup transients and should be discarded.

% Shourov K. Chatterji
% shourov@ligo.caltech.edu
% 2008-Aug-01

% timing parameters
runStart = 894376679;
runDuration = 1024 * 1;%106;
channelName = 'H1:STRAIN';
frameType = 'H1_NINJA_NOISE';
detectorName = channelName(1:2);
siteName = channelName(1);
frameStride = 1024;
sampleFrequency = 4096;
nyquistFrequency = sampleFrequency / 2;
timeResolution = 1 / sampleFrequency;
time = 0 : timeResolution : frameStride - timeResolution;

% low frequency seismic noise shaping filter
seismicSampleFactor = 1024 / nyquistFrequency;
[zSeismic, pSeismic, kSeismic] = ...
    cheby1(8, 20, 20 / (seismicSampleFactor * nyquistFrequency), 'low');
kSeismic = 1 * 5e-31;
[bSeismic, aSeismic] = zp2tf(zSeismic, pSeismic, kSeismic);

% thermal noise shaping filter
thermalCutoffFrequency = 25;
zThermal = [-1; -1;];
pThermal = [+1; +1;] * (1 -  thermalCutoffFrequency / (nyquistFrequency / 2));
kThermal = 1 * 1.0e-23;
[bThermal, aThermal] = zp2tf(zThermal, pThermal, kThermal);

% shot noise shaping filter
shotSampleFactor = 2;
shotCutoffFrequency = 1 * 200;
zShot = [+1; +1 - shotCutoffFrequency / (nyquistFrequency * shotSampleFactor / 2)];
pShot = [-1; +1];
kShot = 1 * 2.0e-20;
[bShot, aShot] = zp2tf(zShot, pShot, kShot);

% anti-alias filter
[zAntiAlias, pAntiAlias, kAntiAlias] = ...
    butter(20, 0.97 / shotSampleFactor, 'low');
kAntiAlias = kAntiAlias;
[bAntiAlias, aAntiAlias] = zp2tf(zAntiAlias, pAntiAlias, kAntiAlias);

% number of frames to produce
numberOfFrames = ceil(runDuration / frameStride);

% seismic noise interpolation filter
seismicInterpolater = intfilt(1 / seismicSampleFactor, ...
                              4 / seismicSampleFactor, ...
                              0.9);

% initialize filter state vectors
zSeismic = [];
zThermal = [];
zShot = [];
zAntiAlias = [];
zSeismicInterp = [];

% gps start times of frames
frameTimes = runStart + (-1 : numberOfFrames - 1) * frameStride;

% create output directory
[ignore, ignore] = mkdir(['frames/noise/' detectorName]);

% initialize random number generator
randn('state', 100000);

% begin loop over frames
for frameTime = frameTimes,

  % display status
  fprintf(1, 'creating %s-%s-%d-%d.gwf...\n', siteName, frameType, ...
          round(frameTime), round(frameStride));

  % stationary gaussian white noise excitation
  seismicExcitation = randn(1, seismicSampleFactor * sampleFrequency * ...
                            frameStride);
  thermalExcitation = randn(1, sampleFrequency * ...
                            frameStride);
  shotExcitation    = randn(1, shotSampleFactor * sampleFrequency * ...
                            frameStride);

  % filter excitation noise with shaping filters
  [seismicNoise, zSeismic] = filter(bSeismic, aSeismic, ...
                                    seismicExcitation, zSeismic);
  [thermalNoise, zThermal] = filter(bThermal, aThermal, ...
                                    thermalExcitation, zThermal);
  [shotNoise, zShot]       = filter(bShot,    aShot,    ...
                                    shotExcitation, zShot);

  % filter shot noise with with anti-alias filter
  [shotNoise, zAntiAlias] = filter(bAntiAlias, aAntiAlias, shotNoise, ...
                                   zAntiAlias);

  % downsample shot noise to target sample frequency
  shotNoise = shotNoise(1 : shotSampleFactor : end);

  % upsample seismic noise to target sample frequency
  seismicNoise = upsample(seismicNoise, 1 / seismicSampleFactor);
  [seismicNoise, zSeismicInterp] = filter(seismicInterpolater, 1, ...
                                          seismicNoise, zSeismicInterp);

  % total detector noise
  detectorNoise = seismicNoise + thermalNoise + shotNoise;

  % prepare structure for frame file
  frameData.channel = channelName;
  frameData.data = detectorNoise;
  frameData.type = 'd';
  frameData.mode = 'a';

  % frame file name
  frameFile = sprintf('frames/noise/%s/%s-%s-%d-%d.gwf', ...
                      detectorName, siteName, frameType, ...
                      round(frameTime), round(frameStride));

  % write frame file
  if frameTime >= runStart,
    mkframe(frameFile, frameData, 'n', frameStride, frameTime);
  end

% end loop over frames
end

% report status
disp('done.');

% return to calling function
return;
