modulated.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2016 - 2021 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
// epsic/src/modulated.h
 
#ifndef __epsic_modulated_h
#define __epsic_modulated_h
 
#include "mode.h"
 
#include <vector>
 
#define _DEBUG 0
#if _DEBUG
#include <iostream>
#endif
 
namespace epsic
{
 
  /***************************************************************************
   *
   *  an amplitude modulated mode of electromagnetic radiation
   *
   ***************************************************************************/
 
  class modulated_mode : public field_transformer
  {
  #if _DEBUG
    mutable double tot, totsq;
    uint64_t count;
  #endif
 
  public:
 
    modulated_mode (mode* s) : field_transformer(s)
    { 
  #if -DEBUG
      tot=0; totsq=0; count=0; 
  #endif
    }
 
    // return a random scalar modulation factor
    virtual double modulation () = 0;
 
    // return the mean of the scalar modulation factor
    virtual double get_mod_mean () const = 0;
 
    // return the variance of the scalar modulation factor
    virtual double get_mod_variance () const = 0;
 
    Spinor<double> transform (const Spinor<double>& field)
    {
      double mod = modulation();
  #if _DEBUG
      tot+=mod;
      totsq+=mod*mod;
      count+=1;
  #endif
      return sqrt(mod) * field;
    }
 
    Matrix<4,4,double> get_covariance () const
    {
      double mean = get_mod_mean();
      double var = get_mod_variance();
 
  #if _DEBUG
      cerr << "modulated_mode::get_covariance expected mean=" << mean
    << " var=" << var << endl;
      tot /= count;
      totsq /= count;
      totsq -= tot*tot;
      cerr << " measured mean=" << tot << " var=" << totsq << endl;
  #endif
 
      Matrix<4,4,double> C = source->get_covariance();
      C *= (mean*mean + var);
      Matrix<4,4,double> o = outer (source->get_mean(), source->get_mean());
      o *= var;
      return C + o;
    }
 
    Stokes<double> get_mean () const
    {
      return get_mod_mean () * source->get_mean();
    }
 
  };
 
  class lognormal_mode : public modulated_mode
  {
    // standard deviation of the logarithm of the random variate
    double log_sigma;
 
  public:
 
    // initialize based on the modulation index beta
    lognormal_mode (mode* s, double beta) : modulated_mode (s) { set_beta (beta); }
 
    void set_beta (double beta)
    {
      log_sigma = sqrt( log( beta*beta + 1.0 ) );
    }
 
    double get_beta () const
    {
      return sqrt( get_mod_variance() );
    }

    double get_log_sigma () const
    {
      return log_sigma;
    }
 
    // return a random scalar modulation factor
    double modulation ()
    {
      return exp ( log_sigma * (get_normal()->evaluate() - 0.5*log_sigma) ) ;
    }
 
    double get_mod_mean () const { return 1.0; }
    
    double get_mod_variance () const { return exp(log_sigma*log_sigma) - 1.0; }
 
  };
 
 
  class boxcar_modulated_mode : public modulated_mode
  {
    std::vector< double > instances;
    unsigned smooth;
    unsigned current;
 
    void setup()
    {
      current = 0;
      instances.resize (smooth);
      for (unsigned i=1; i<smooth; i++)
        instances[i] = mod->modulation();
    }
 
    modulated_mode* mod;
 
  public:
 
    boxcar_modulated_mode (modulated_mode* s, unsigned n)
      : modulated_mode(s->get_source()) { smooth = n; mod = s; }
 
    double modulation ()
    {
      if (instances.size() < smooth)
        setup ();
 
      instances[current] = mod->modulation();
      current = (current + 1) % smooth;
 
      double result = 0.0;
      for (unsigned i=0; i<smooth; i++)
        result += instances[i];
 
      result /= smooth;
 
      return result;
    }
 
    double get_mod_variance () const
    {
      return mod->get_mod_variance() / smooth;
    }
 
    double get_mod_mean () const
    {
      return mod->get_mod_mean ();
    }
 
    Matrix<4,4, double> get_crosscovariance (unsigned ilag) const
    {
      if (ilag >= smooth)
        return 0;
 
      Matrix<4,4, double> result = outer(source->get_mean(), source->get_mean());
      result *= double (smooth - ilag) / smooth * get_mod_variance();
      return result;
    }
    
  };
 
 
  class square_modulated_mode : public modulated_mode
  {
    unsigned width;
    unsigned current;
    double value;
 
    modulated_mode* mod;
 
    std::vector<double> cross_correlation;
 
  public:
 
    square_modulated_mode (modulated_mode*,
        unsigned width,
        unsigned sample_size);
 
    void compute_cross_correlation (unsigned sample_size);
    
    double modulation ()
    {
      if (current == width)
      {
        value = mod->modulation();
        current = 0;
      }
 
      current ++;
      return value;
    }
 
    double get_mod_variance () const
    {
      return mod->get_mod_variance();
    }
 
    double get_mod_mean () const
    {
      return mod->get_mod_mean ();
    }

    Matrix<4,4, double> get_crosscovariance (unsigned ilag) const
    {
      if (ilag >= width)
        return 0;
 
      Matrix<4,4, double> result = outer(source->get_mean(), source->get_mean());
      result *= cross_correlation[ilag] * get_mod_variance();
 
      return result;
    }
  };
 
} // namespace epsic
 
#endif // ! defined __epsic_modulated_h