distmain.cpp

Digital filter designer s handbook C++ code source

//
//  File = distmain.cpp
//
//
//  Program for Probability Distributions
//
#include <fstream.h>
#include <math.h>  
#include "misdefs.h"
#include "uni_rand.h"
#include "gausrand.h"
#include "gaus_pdf.h"
#include "exp_rand.h"
#include "raylrand.h"
#include "srand.h"

ofstream DebugFile("distrib.bug", ios::out);
ofstream PdfFile("gaus_pdf.txt", ios::out); 
ofstream CdfFile("gaus_cdf.txt", ios::out); 
ofstream HistoFile("gaushist.txt", ios::out);

main()
{
 long seed;
 long num_samps;
 long samp_indx;
 long *bin_counts;
 int num_bins, bin_indx, i;
 float ran_samp;
 double double_ran_samp;
 int half_num_bins;
 int num_bins_div_8;
 double g1;
 double std_dev;
 int int_g1, int_seed;
 int plot_type, distrib_type, num_plot_pts;
 srand* PNG_1;
 ofstream dec_file("dec_exp.txt", ios::out); 
 ofstream sat_file("sat_exp.txt", ios::out); 

 cout << "Desired type of plot?\n"
      << "  1 = empirical PDF\n"
      << "  2 = empirical CDF\n"
      << "  3 = theoretical PDF\n"
      << "  4 = theoretical CDF\n"
      << "  5 = decaying exponential\n"
      << "  6 = saturating exponential\n"
      << endl;
 cin >> plot_type;
 
 switch(plot_type)
  {
  case 1:   //empirical PDF
  case 2:   //empirical CDF
    cout << "Desired distribution?\n"
         << "  1 = uniform\n"
         << "  2 = gaussian\n"
         << "  3 = exponential\n"
         << "  4 = rayleigh\n"
         << endl;
    cin >> distrib_type;
    cout << "starting seed?" << endl;
    cin >> seed;
    cout << "how many samples?" << endl;
    cin >> num_samps;
    cout << "how many bins in histogram?" << endl;
    cin >> num_bins;
    half_num_bins = num_bins/2;
    num_bins_div_8 = num_bins/8;
    bin_counts = new long[num_bins];
    for(i=0; i<num_bins; i++) bin_counts[i] = 0; 
    
    int_seed = (unsigned int)(seed & 32767);
    PNG_1 = new srand(int_seed);
 
    for(samp_indx=0; samp_indx<num_samps; samp_indx++)
      {
       switch (distrib_type)
        {
        case 1:    //uniform
          //g1 = DoubleUniformRandom(&seed);
          int_g1 = PNG_1->rand();
          //cout << int_g1 << endl;
          g1 = float(int_g1)/float(RAND_MAX);
          //cout << g1 << endl;
          break;
        case 2:    // gaussian
          g1 = GaussRandom(&seed);
          break;
        case 3:    // exponential
          g1 = ExponentialRandom(&seed);
          break;
        case 4:    // Rayleigh
          g1 = RayleighRandom(&seed);
          break;
        }

       //bin_indx = half_num_bins + int( num_bins_div_8 * g1 );
       bin_indx = int( num_bins * g1 );
       if(g1 < 0.0) bin_indx--;
       if(bin_indx < 0) continue;
       if(bin_indx >= num_bins) continue;
       bin_counts[bin_indx]++;
      }
    for(bin_indx=0; bin_indx<num_bins; bin_indx++)
      {
       HistoFile << bin_indx << ", " 
                 << float(bin_counts[bin_indx])  //float(num_samps)
                 << endl;
      }
    break;
  //------------------------------------------------------   
  case 3:   // theoretical PDF
    cout << "Desired distribution?\n"
         << "  1 = uniform\n"
         << "  2 = gaussian\n"
         << "  3 = exponential\n"
         << "  4 = rayleigh\n"
         << endl;
    cin >> distrib_type;
    cout << "how many points in plot?" << endl;
    cin >> num_plot_pts;
    cout << "desired value of standard deviation?" << endl;
    cin >> std_dev;
    switch (distrib_type)
      {
      case 1:    //uniform
        //UniformPdf();
        break;
      case 2:    // gaussian
        GaussianPdf(std_dev, 0.0, num_plot_pts, 4.0);
        break;
      case 3:    // exponential
        //ExponentialPdf();
        break;
      case 4:    // Rayleigh
        //RayleighPdf();
        break;
      }
    break;
  case 4:  // theoretical CDF
    cout << "Desired distribution?\n"
         << "  1 = uniform\n"
         << "  2 = gaussian\n"
         << "  3 = exponential\n"
         << "  4 = rayleigh\n"
         << endl;
    cin >> distrib_type;
    cout << "how many points in plot?" << endl;
    cin >> num_plot_pts;

    switch (distrib_type)
      {
      case 1:    //uniform
        //UniformCdf();
        break;
      case 2:    // gaussian
        GaussianCdf(1.0, 0.0, num_plot_pts, 4.0);
        break;
      case 3:    // exponential
        //ExponentialCdf();
        break;
      case 4:    // Rayleigh
        //RayleighCdf();
        break;
      }
    break;
  case 5:   // decaying exponential
    double alpha, beta, delta_time, time, y_val;
    int plot_indx;
    cout << "desired value for beta?" << endl;
    cin >> beta;
    cout << "desired value for alpha?" << endl;
    cin >> alpha;
    cout << "number of points in plot?" << endl;
    cin >> num_plot_pts;
    delta_time = 5.0/(alpha * num_plot_pts);
    for(plot_indx=0; plot_indx<num_plot_pts; plot_indx++)
      {
      time = delta_time * plot_indx;
      y_val = beta * exp(-alpha * time);
      dec_file << time << ", " << y_val << endl;
      }
    break;
  case 6:   // saturating exponential
    cout << "desired value for beta?" << endl;
    cin >> beta;
    cout << "desired value for alpha?" << endl;
    cin >> alpha;
    cout << "number of points in plot?" << endl;
    cin >> num_plot_pts;
    delta_time = 5.0/(alpha * num_plot_pts);
    for(plot_indx=0; plot_indx<num_plot_pts; plot_indx++)
      {
      time = delta_time * plot_indx;
      y_val = 1.0 - beta * exp(-alpha * time);
      sat_file << time << ", " << y_val << endl;
      }
    break;
  }
 
 return 0;
}