src/s_fo.cc

   1 /*                             -*- C++ -*-
   2  * Copyright (C) 2001 Albert Davis
   3  * Author: Albert Davis <aldavis@gnu.org>
   4  *
   5  * This file is part of "Gnucap", the Gnu Circuit Analysis Package
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  20  * 02110-1301, USA.
  21  *------------------------------------------------------------------
  22  * tran and fourier commands -- top
  23  * performs transient analysis, silently, then fft.
  24  * outputs results of fft
  25  */
  26 #include "globals.h"
  27 #include "u_sim_data.h"
  28 #include "u_status.h"
  29 #include "m_phase.h"
  30 #include "declare.h"    /* plclose, plclear, fft */
  31 #include "u_prblst.h"
  32 #include "s_tr.h"
  33 /*--------------------------------------------------------------------------*/
  34 namespace {
  35 /*--------------------------------------------------------------------------*/
  36 class FOURIER : public TRANSIENT {
  37 public:
  38   void  do_it(CS&, CARD_LIST*);
  39   explicit FOURIER():
  40     TRANSIENT(),
  41     _fstart(0.),
  42     _fstop(0.),
  43     _fstep(0.),
  44     _timesteps(0),
  45     _fdata(NULL)
  46   {}
  47   ~FOURIER() {}
  48 private:
  49   explicit FOURIER(const FOURIER&): TRANSIENT() {unreachable(); incomplete();}
  50   std::string status()const {untested();return "";}
  51   void  setup(CS&);     /* s_fo_set.cc */
  52   void  fftallocate();
  53   void  fftunallocate();
  54   void  foout();        /* s_fo_out.cc */
  55   void  fohead(const PROBE&);
  56   void  foprint(COMPLEX*);
  57   void  store_results(double); // override virtual
  58 private:
  59   PARAMETER<double> _fstart;    /* user start frequency */
  60   PARAMETER<double> _fstop;     /* user stop frequency */
  61   PARAMETER<double> _fstep;     /* fft frequecncy step */
  62   int    _timesteps;    /* number of time steps in tran analysis, incl 0 */
  63   COMPLEX** _fdata;     /* storage to allow postprocessing */
  64 };
  65 /*--------------------------------------------------------------------------*/
  66 static  int     to_pow_of_2(double);
  67 static  int     stepnum(double,double,double);
  68 static  COMPLEX find_max(COMPLEX*,int,int);
  69 static  double  db(COMPLEX);
  70 /*--------------------------------------------------------------------------*/
  71 void FOURIER::do_it(CS& Cmd, CARD_LIST* Scope)
  72 {
  73   _scope = Scope;
  74   _sim->set_command_fourier();
  75   reset_timers();
  76   ::status.four.reset().start();
  77
  78   try {
  79     setup(Cmd);
  80     _sim->init();
  81     CARD_LIST::card_list.precalc_last();
  82
  83     _sim->alloc_vectors();
  84     _sim->_aa.reallocate();
  85     _sim->_aa.dezero(OPT::gmin);
  86     _sim->_aa.set_min_pivot(OPT::pivtol);
  87     _sim->_lu.reallocate();
  88     _sim->_lu.dezero(OPT::gmin);
  89     _sim->_lu.set_min_pivot(OPT::pivtol);
  90     fftallocate();
  91     ::status.set_up.stop();
  92
  93     switch (ENV::run_mode) {
  94     case rPRE_MAIN:     unreachable();          break;
  95          case rPIPE:
  96     case rBATCH:        untested();
  97     case rINTERACTIVE:  itested();
  98     case rSCRIPT:       sweep(); foout();       break;
  99     case rPRESET:       untested(); /*nothing*/ break;
 100     }
 101   }catch (Exception& e) {untested();
 102     error(bDANGER, e.message() + '\n');
 103   }
 104   fftunallocate();
 105   _sim->unalloc_vectors();
 106   _sim->_lu.unallocate();
 107   _sim->_aa.unallocate();
 108
 109   _sim->_has_op = s_FOURIER;
 110   _scope = NULL;
 111
 112   ::status.four.stop();
 113   ::status.total.stop();
 114
 115 }
 116 /*--------------------------------------------------------------------------*/
 117 /*--------------------------------------------------------------------------*/
 118 /* store: stash time domain data in preparation for Fourier Transform
 119  */
 120 void FOURIER::store_results(double X)
 121 {
 122   TRANSIENT::store_results(X);
 123
 124   if (step_cause() == scUSER) {
 125     int ii = 0;
 126     for (PROBELIST::const_iterator
 127            p=printlist().begin();  p!=printlist().end();  ++p) {
 128       assert(_sim->_stepno < _timesteps);
 129       _fdata[ii][_sim->_stepno] = (*p)->value();
 130       ++ii;
 131     }
 132   }else{untested();
 133   }
 134 }
 135 /*--------------------------------------------------------------------------*/
 136 /* foout:  print out the results of the transform
 137  */
 138 void FOURIER::foout()
 139 {
 140   plclose();
 141   plclear();
 142   int ii = 0;
 143   for (PROBELIST::const_iterator
 144          p=printlist().begin();  p!=printlist().end();  ++p) {
 145     fohead(**p);
 146     fft(_fdata[ii], _timesteps-1,  0);
 147     foprint(_fdata[ii]);
 148     ++ii;
 149   }
 150 }
 151 /*--------------------------------------------------------------------------*/
 152 /* fo_head: print output header
 153  * arg is index into probe array, to select probe name
 154  */
 155 void FOURIER::fohead(const PROBE& Prob)
 156 {
 157   _out.form("# %-10s", Prob.label().c_str())
 158     << "--------- actual ---------  -------- relative --------\n"
 159     << "#freq       "
 160     << "value        dB      phase  value        dB      phase\n";
 161 }
 162 /*--------------------------------------------------------------------------*/
 163 /* fo_print: print results of fourier analysis
 164  * for all points at single probe
 165  */
 166 void FOURIER::foprint(COMPLEX *Data)
 167 {
 168   int startstep = stepnum(0., _fstep, _fstart);
 169   assert(startstep >= 0);
 170   int stopstep  = stepnum(0., _fstep, _fstop );
 171   assert(stopstep < _timesteps);
 172   COMPLEX maxvalue = find_max(Data, std::max(1,startstep), stopstep);
 173   if (maxvalue == 0.) {untested();
 174     maxvalue = 1.;
 175   }else{
 176   }
 177   Data[0] /= 2;
 178   for (int ii = startstep;  ii <= stopstep;  ++ii) {
 179     double frequency = _fstep * ii;
 180     assert(ii >= 0);
 181     assert(ii < _timesteps);
 182     COMPLEX unscaled = Data[ii];
 183     COMPLEX scaled = unscaled / maxvalue;
 184     unscaled *= 2;
 185     _out.form("%s%s%7.2f %8.3f %s%7.2f %8.3f\n",
 186              ftos(frequency,    11,5,_out.format()),
 187         ftos(std::abs(unscaled),11,5,_out.format()),
 188              db(unscaled),
 189              phase(unscaled*COMPLEX(0.,1)),
 190         ftos(std::abs(scaled),  11,5,_out.format()),
 191              db(scaled),
 192              phase(scaled) ) ;
 193   }
 194 }
 195 /*--------------------------------------------------------------------------*/
 196 /* stepnum: return step number given its frequency or time
 197  */
 198 static int stepnum(double Start, double Step, double Here)
 199 {
 200   return int((Here-Start)/Step + .5);
 201 }
 202 /*--------------------------------------------------------------------------*/
 203 /* find_max: find the max magnitude in a COMPLEX array
 204  */
 205 static COMPLEX find_max(COMPLEX *Data, int Start, int Stop)
 206 {
 207   COMPLEX maxvalue = 0.;
 208   for (int ii = Start;  ii <= Stop;  ++ii) {
 209     if (std::abs(Data[ii]) > std::abs(maxvalue)) {
 210       maxvalue = Data[ii];
 211     }else{
 212     }
 213   }
 214   return maxvalue;
 215 }
 216 /*--------------------------------------------------------------------------*/
 217 static double db(COMPLEX Value)
 218 {
 219   return  20. * log10(std::max(std::abs(Value),VOLTMIN));
 220 }
 221 /*--------------------------------------------------------------------------*/
 222 /*--------------------------------------------------------------------------*/
 223 /* fo_setup: fourier analysis: parse command string and set options
 224  *      (options set by call to TRANSIENT::options)
 225  */
 226 void FOURIER::setup(CS& Cmd)
 227 {
 228   _cont = true;
 229   if (Cmd.match1("'\"({") || Cmd.is_pfloat()) {
 230     PARAMETER<double> arg1, arg2, arg3;
 231     Cmd >> arg1;
 232     if (Cmd.match1("'\"({") || Cmd.is_float()) {
 233       Cmd >> arg2;
 234     }else{untested();
 235     }
 236     if (Cmd.match1("'\"({") || Cmd.is_float()) {
 237       Cmd >> arg3;
 238     }else{untested();
 239     }
 240
 241     if (arg3.has_hard_value()) {            /* 3 args: all */
 242       assert(arg2.has_hard_value());
 243       assert(arg1.has_hard_value());
 244       _fstart = arg1;
 245       _fstop  = arg2;
 246       _fstep  = arg3;
 247     }else if (arg2.has_hard_value()) {untested(); /* 2 args: start = 0 */
 248       assert(arg1.has_hard_value());
 249       arg1.e_val(0.,_scope);
 250       arg2.e_val(0.,_scope);
 251       if (arg1 >= arg2) {untested();        /* 2 args: stop, step */
 252         _fstart = "NA";                     /*  (stop > step) */
 253         _fstop  = arg1;
 254         _fstep  = arg2;
 255       }else{untested(); /* arg1 < arg2 */    /* 2 args: step, stop */
 256         _fstart = "NA";
 257         _fstop  = arg2;
 258         _fstep  = arg1;
 259       }
 260     }else{untested();
 261       assert(arg1.has_hard_value());
 262       arg1.e_val(0.,_scope);
 263       if (arg1 == 0.) {untested();          /* 1 arg: start */
 264         _fstart = 0.;
 265         /* _fstop unchanged */
 266         /* _fstep unchanged */
 267       }else{untested();                     /* 1 arg: step */
 268         _fstart = "NA";
 269         _fstop  = "NA";
 270         _fstep  = arg1;
 271       }
 272     }
 273   }else{untested();
 274     /* else (no args) : no change */
 275   }
 276
 277   options(Cmd);
 278
 279   _fstart.e_val(0., _scope);
 280   _fstep.e_val(0., _scope);
 281   _fstop.e_val(OPT::harmonics * _fstep, _scope);
 282
 283   if (_fstep == 0.) {untested();
 284     throw Exception("frequency step = 0");
 285   }else{
 286   }
 287   if (_fstop == 0.) {untested();
 288     _fstop = OPT::harmonics * _fstep;
 289   }else{
 290   }
 291
 292   _timesteps = to_pow_of_2(_fstop*2 / _fstep) + 1;
 293   if (_cold  ||  _sim->last_time() <= 0.) {
 294     _cont = false;
 295     _tstart = 0.;
 296   }else{
 297     _cont = true;
 298     _tstart = _sim->last_time();
 299   }
 300   _tstop = _tstart + 1. / _fstep;
 301   _tstrobe = 1. / _fstep / (_timesteps-1);
 302   _time1 = _sim->_time0 = _tstart;
 303
 304   _sim->_freq = _fstep;
 305
 306   _dtmax = std::min(double(_dtmax_in), _tstrobe / double(_skip_in));
 307   if (_dtmin_in.has_hard_value()) {untested();
 308     _sim->_dtmin = _dtmin_in;
 309   }else if (_dtratio_in.has_hard_value()) {untested();
 310     _sim->_dtmin = _dtmax / _dtratio_in;
 311   }else{
 312     // use smaller of soft values
 313     _sim->_dtmin = std::min(double(_dtmin_in), _dtmax/_dtratio_in);
 314   }
 315 }
 316 /*--------------------------------------------------------------------------*/
 317 /* allocate:  allocate space for fft
 318  */
 319 void FOURIER::fftallocate()
 320 {
 321   int probes = printlist().size();
 322   _fdata = new COMPLEX*[probes];
 323   for (int ii = 0;  ii < probes;  ++ii) {
 324     _fdata[ii] = new COMPLEX[_timesteps+100];
 325   }
 326 }
 327 /*--------------------------------------------------------------------------*/
 328 /* unallocate:  unallocate space for fft
 329  */
 330 void FOURIER::fftunallocate()
 331 {
 332   if (_fdata) {
 333     for (int ii = 0;  ii < printlist().size();  ++ii) {
 334       delete [] _fdata[ii];
 335     }
 336     delete [] _fdata;
 337     _fdata = NULL;
 338   }else{untested();
 339   }
 340 }
 341 /*--------------------------------------------------------------------------*/
 342 /* to_pow_of_2: round up to nearest power of 2
 343  * example: z=92 returns 128
 344  */
 345 static int to_pow_of_2(double Z)
 346 {
 347   int x = static_cast<int>(floor(Z));
 348   int y;
 349   for (y = 1; x > 0; x >>= 1) {
 350     y <<= 1;
 351   }
 352   return y;
 353 }
 354 /*--------------------------------------------------------------------------*/
 355 static FOURIER p3;
 356 DISPATCHER<CMD>::INSTALL d3(&command_dispatcher, "fourier", &p3);
 357 }
 358 /*--------------------------------------------------------------------------*/
 359 /*--------------------------------------------------------------------------*/
 360 // vim:ts=8:sw=2:noet: