// $Id$
//
// On-line normalization.  This block is intended to replace our
// standard off-line normalization in which means and variances are
// calculated over an entire data set and then applied to the data
// prior to MLP training to generate zero-mean, unit-variance
// features.  This block uses a single-pole low-pass filter
// (with user-specifiable time constant) to estimate the signal mean
// on-line.  This mean is subtracted off, then a second single
// low-pass filter (again, with user-specifiable time constant),
// preceded by a squaring operation, is used to estimate the signal
// variance.  The zero-mean signal is divided by the square root of
// this estimate, plus a small epsilon to prevent divide-by-zero
// problems.  Because the filters used have relatively long time
// constants and we may be processing relatively short files, a
// history mechanism is employed so that the filter state from a
// previous utterance may be used to initialize the filters for the
// current utterance.  On termination, the filter state may be written
// to a history file for use by later runs.
//
// $Log$
//

#include <assert.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include "Online_Norm.h"

const CMatrix 
Online_Norm::processBlock(const CMatrix in,
                          bool& lastTime)
{
  assert(in.m == numChannels);
  assert(in.n <= IOBlocking);

  if (state == End)
    {
      lastTime = true;
      return CMatrix(numChannels,
                     0,
                     IOBlocking,
                     result);
    }

  if (in.n == 0)
    {
      if (lastTime) { state = End; }
      return CMatrix(numChannels,
                     0,
                     IOBlocking,
                     result);
    }

  // variables for accessing input matrix
  //
  const size_t in_row = in.m;
  const size_t in_col = in.n;
  const size_t in_stride = in.s;
  const float* in_row_startp = in.data;
  const float* in_row_endp = in.data + in_col;
  const float* in_endp = in.data + in_row * in_stride;

  // variables for accessing result.
  //
  float* res_row_startp = result;

  float* mean_p = meanLPstate;
  float* var_p = varLPstate;

  while (in_row_startp != in_endp)
    {
      const float* in_p = in_row_startp;
      float* res_p = res_row_startp;
      float y = *mean_p;
      float z = *var_p;
      while (in_p != in_row_endp)
        {
          float x = *in_p++;
          y = meanLP_a1 * y + meanLP_b0 * x;
          x = x - y;
          z = varLP_a1 * z + varLP_b0 * x * x;
          x = x / (sqrt(z) + epsilon);
          *res_p++ = x;
        }
      *mean_p++ = y;
      *var_p++ = z;
      in_row_startp += in_stride;
      in_row_endp += in_stride;
      res_row_startp += IOBlocking;
    }

  if (lastTime) { state = End; }

  return CMatrix(numChannels,
                 in_col,
                 IOBlocking,
                 result);
}

void 
Online_Norm::readHistory()
{
  size_t nc;
  FILE* inHist;
  char magic[8];

  if ((inHist = fopen(inHistoryFileName, "r"))  == NULL)
    {
      fprintf(stderr,
              "Unable to open %s for input: aborting.\n",
              inHistoryFileName);
      exit(EXIT_FAILURE);
    }

  if (fread(magic, sizeof(char), 7, inHist) != 7)
    {
      fprintf(stderr,
              "Error reading %s: aborting.\n",
              inHistoryFileName);
      exit(EXIT_FAILURE);
    }
  magic[7] = '\000';
  if (strcmp(magic, "HISTORY") != 0)
    {
      fprintf(stderr,
              "%s is not a valid history file: aborting.\n",
              inHistoryFileName);
      exit(EXIT_FAILURE);
    }

  if (fread(&nc, sizeof(size_t), 1, inHist) != 1)
    {
      fprintf(stderr,
              "Error reading %s: aborting.\n",
              inHistoryFileName);
      exit(EXIT_FAILURE);
    }
  if (nc != numChannels)
    {
      fprintf(stderr,
              "%s is for processing %d channels, but we are processing %d channels: aborting.\n",
              inHistoryFileName,
              nc,
              numChannels);
      exit(EXIT_FAILURE);
    }

  if (fread(meanLPstate, sizeof(float), numChannels, inHist) != numChannels)
    {
      fprintf(stderr,
              "Error reading %s: aborting.\n",
              inHistoryFileName);
      exit(EXIT_FAILURE);
    }

  if (fread(varLPstate, sizeof(float), numChannels, inHist) != numChannels)
    {
      fprintf(stderr,
              "Error reading %s: aborting.\n",
              inHistoryFileName);
      exit(EXIT_FAILURE);
    }

  if (fclose(inHist) == EOF)
    {
      fprintf(stderr,
              "Error closing %s: aborting.\n",
              inHistoryFileName);
      exit(EXIT_FAILURE);
    }
}

void 
Online_Norm::writeHistory()
{
  FILE* outHist;
  const char* magic = "HISTORY";

  if ((outHist = fopen(outHistoryFileName, "w"))  == NULL)
    {
      fprintf(stderr,
              "Unable to open %s for output: aborting.\n",
              outHistoryFileName);
      exit(EXIT_FAILURE);
    }

  if (fwrite(magic, sizeof(char), 7, outHist) != 7)
    {
      fprintf(stderr,
              "Error writing %s: aborting.\n",
              outHistoryFileName);
      exit(EXIT_FAILURE);
    }

  if (fwrite(&numChannels, sizeof(size_t), 1, outHist) != 1)
    {
      fprintf(stderr,
              "Error writing %s: aborting.\n",
              outHistoryFileName);
      exit(EXIT_FAILURE);
    }

  if (fwrite(meanLPstate, sizeof(float), numChannels, outHist) != numChannels)
    {
      fprintf(stderr,
              "Error writing %s: aborting.\n",
              outHistoryFileName);
      exit(EXIT_FAILURE);
    }

  if (fwrite(varLPstate, sizeof(float), numChannels, outHist) != numChannels)
    {
      fprintf(stderr,
              "Error writing %s: aborting.\n",
              outHistoryFileName);
      exit(EXIT_FAILURE);
    }

  if (fclose(outHist) == EOF)
    {
      fprintf(stderr,
              "Error closing %s: aborting.\n",
              outHistoryFileName);
      exit(EXIT_FAILURE);
    }
}