home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c100 / 2.ddi / SEQUEN.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1988-08-20  |  11.0 KB  |  397 lines
  1. #include <iostream.h>
  2. #include <alloc.h>
  3. #include <stdlib.h>
  4.  
  5. /*---------------04-26-90 09:29am-------------------
  6.  * An illustration of defining and using a class which employs
  7.  * dynamically allocated arrays.  The sequence class is used to manipulate
  8.  * time-domain sequences for digital signal processing.  We have an array
  9.  * of signal values, and the time extent (in samples) which the signal
  10.  * encompasses.  We define operators for storing (=), convolving (*), adding,
  11.  * subtracting, IIR convolution (/), construction and destruction.
  12.  --------------------------------------------------*/
  13. class sequence {
  14.     protected:
  15.  
  16.         double huge *seq;  //sequence of points - some signals are long
  17.         int lower_limit;   //start time
  18.         int upper_limit;   //end time
  19.  
  20.     public:
  21.         sequence(int lower=0, int upper=0);  //default sequence is single point
  22.         sequence(sequence& arg);  //construct from other sequence
  23.         ~sequence(void);  //destructor
  24.         double fetch(int index);  //need operator to get datum
  25.         void stuff(int index,double val);  //and an operator to place datum
  26.         void extent(int *lower, int *upper); //since limits are hidden, fetch them
  27.         sequence operator* (sequence& arg); //convolution
  28.         sequence& operator*= (sequence& arg);  //convolution
  29.         sequence operator/ (sequence& arg); //IIR filter
  30.         sequence& operator/= (sequence& arg);  //IIR filter
  31.         sequence operator+ (sequence& arg);  //add sequences
  32.         sequence& operator+= (sequence &arg); //add in place
  33.         sequence operator- (sequence& arg);  //subtract
  34.         sequence& operator-= (sequence &arg);  //subtract in place
  35.         sequence& operator= (sequence& arg);  //copy
  36.         int is_causal(void);  //necessary for DSP
  37. };
  38.  
  39. /*create a sequence of zeros, from lower to upper index*/
  40. sequence::sequence(int lower, int upper)
  41. {
  42.     int i;
  43.  
  44.     if (lower>upper) {
  45.         i=lower;
  46.         lower=upper;
  47.         upper=i;
  48.     }
  49.     seq=new double[upper+1-lower];
  50.     if (seq==NULL) {
  51.         cerr << "Cannot allocate sequence\n";
  52.         exit(1);
  53.     }
  54.     lower_limit=lower;
  55.     upper_limit=upper;
  56.     for (i=lower;i<=upper;i++) {
  57.         seq[i-lower]=0.0;
  58.     }
  59. }
  60.  
  61. /*create a sequence from another sequence*/
  62. sequence::sequence(sequence& arg)
  63. {
  64.     int i;
  65.  
  66.     /*set limits*/
  67.     lower_limit=arg.lower_limit;
  68.     upper_limit=arg.upper_limit;
  69.     /*allocate the array*/
  70.     seq=new double[upper_limit+1-lower_limit];
  71.     if(seq==NULL){
  72.         cerr << "Cannot allocate sequence\n";
  73.         exit(1);
  74.     }
  75.     /*transfer the data*/
  76.     for(i=lower_limit;i<=upper_limit;i++){
  77.         stuff(i,arg.fetch(i));
  78.     }
  79. }
  80.  
  81. /*destroy a sequence*/
  82. sequence::~sequence(void)
  83. {
  84.     delete seq;
  85. }
  86.  
  87. /*since the data array is hidden, we need this to put data in, if valid*/
  88. void sequence::stuff(int index,double val)
  89. {
  90.     if ((index>=lower_limit)&&(index<=upper_limit)) {
  91.         seq[index-lower_limit]=val;
  92.     }
  93. }
  94.  
  95. /*since the data array is hidden, we need this to put data in, if valid*/
  96. double sequence::fetch(int index)
  97. {
  98.     if ((index>=lower_limit)&&(index<=upper_limit)) {
  99.         return(seq[index-lower_limit]);
  100.     } else {
  101.         return((double) 0);
  102.     }
  103. }
  104.  
  105. /*since the limits are hidden, we use this to obtain them*/
  106. void sequence::extent(int *lower, int *upper)
  107. {
  108.     *lower=lower_limit;
  109.     *upper=upper_limit;
  110. }
  111.  
  112. /************************************************
  113.  *   --- sequence sequence::operator*  ---
  114.  *    How to convolve two sequences.  NOTE that the operator is NOT by
  115.  * reference.  This is necessary, so as to NOT modify the left-hand
  116.  * sequence
  117.  ************************************************/
  118. sequence sequence::operator* (sequence& arg)
  119. {
  120.     int i,j,lower=lower_limit+arg.lower_limit;
  121.     int upper=upper_limit+arg.upper_limit;
  122.     int this_bigger=(upper_limit+arg.lower_limit-
  123.             (lower_limit+arg.upper_limit));
  124.     /*---------------04-26-90 09:41am-------------------
  125.  * The output is by value.  The temp sequence is used to hold
  126.      * the value
  127.      --------------------------------------------------*/
  128.     sequence temp(lower,upper);
  129.  
  130.     /*To speed the operation, we can reverse the order of the
  131.      * Convolution, depending upon which sequence is longer*/
  132.     if(this_bigger>0){
  133.         /*Convoltion sum one direction*/
  134.         for (i=lower;i<=upper;i++) {
  135.             temp.stuff(i,0.0);
  136.             for (j=arg.lower_limit;j<=arg.upper_limit;j++) {
  137.                 temp.stuff(i,temp.fetch(i)+fetch(i-j)*arg.fetch(j));
  138.             }
  139.         }
  140.     } else {
  141.         /*convolution sum other direction*/
  142.         for (i=lower;i<=upper;i++) {
  143.             temp.stuff(i,0.0);
  144.             for (j=lower_limit;j<=upper_limit;j++) {
  145.                 temp.stuff(i,temp.fetch(i)+fetch(j)*arg.fetch(i-j));
  146.             }
  147.         }
  148.     }
  149.     return temp;
  150. }
  151.  
  152. /*in-place convolution (sort of)*/
  153. sequence& sequence::operator*= (sequence& arg)
  154. {
  155.     sequence temp(*this);
  156.  
  157.     *this=temp*arg;
  158.     return *this;
  159. }
  160.  
  161. /************************************************
  162.  *   --- sequence sequence::operator/  ---
  163.  *    
  164.  *    How to IIR filter.  NOTE that the operator is NOT by
  165.  * reference.  This is necessary, so as to NOT modify the left-hand
  166.  * sequence
  167.  ************************************************/
  168. sequence sequence::operator/ (sequence& arg)
  169. {
  170.     int i,j;
  171.     /*the TEMP sequence is used to hold the output for value output*/
  172.     sequence temp(lower_limit,upper_limit);
  173.     double tempd1,tempd2;
  174.  
  175.     /*We can only do IIR filter if right-hand signal is causal*/
  176.     if((!arg.is_causal())||(arg.fetch(0)==0.0)){
  177.         cerr << "Divide by non-causal sequence\n";
  178.         exit(1);
  179.     }
  180.     /*perform the IIR filter operation*/
  181.     tempd1=1.0/arg.fetch(0);
  182.     for (i=lower_limit;i<=upper_limit;i++) {
  183.         tempd2=fetch(i);
  184.         for (j=arg.lower_limit+1;j<=arg.upper_limit;j++) {
  185.             tempd2-=temp.fetch(i-j)*arg.fetch(j);
  186.         }
  187.         temp.stuff(i,tempd2*tempd1);
  188.     }
  189.     return temp;
  190. }
  191. /*In-place IIR filtering (really) - NOTE:  this IS call-by-reference*/
  192. sequence& sequence::operator/= (sequence& arg)
  193. {
  194.     int i,j;
  195.     double tempd1,tempd2;
  196.  
  197.     if((!arg.is_causal())||(arg.fetch(0)==0.0)){
  198.         cerr << "Divide by non-causal sequence\n";
  199.         exit(1);
  200.     }
  201.     tempd1=1.0/arg.fetch(0);
  202.     for (i=lower_limit;i<=upper_limit;i++) {
  203.         tempd2=fetch(i);
  204.         for (j=arg.lower_limit+1;j<=arg.upper_limit;j++) {
  205.             tempd2-=fetch(i-j)*arg.fetch(j);
  206.         }
  207.         stuff(i,tempd2*tempd1);
  208.     }
  209.     return *this;
  210. }
  211.  
  212. /************************************************
  213.  *   --- sequence sequence::operator+  ---
  214.  *    Add two sequences
  215.  ************************************************/
  216. sequence sequence::operator+ (sequence& arg)
  217. {
  218.     int i,lower=min(lower_limit,arg.lower_limit);
  219.     int upper=max(upper_limit,arg.upper_limit);
  220.     sequence temp(lower,upper);
  221.  
  222.     for (i=lower;i<=upper;i++) {
  223.         temp.stuff(i,fetch(i)+arg.fetch(i));
  224.     }
  225.     return temp;
  226. }
  227.  
  228. /************************************************
  229.  *  --- sequence& sequence::operator+=  ---
  230.  *    Add in-place.  NOTE:  Call-by-reference
  231.  ************************************************/
  232. sequence& sequence::operator+= (sequence& arg)
  233. {
  234.     if ((lower_limit==arg.lower_limit)&&(upper_limit==arg.upper_limit)) {
  235.         int i,lim=upper_limit-lower_limit;
  236.         for (i=0;i<=lim;i++)seq[i]+=arg.seq[i];
  237.     } else {
  238.         sequence temp(*this);
  239.         *this=temp+arg;
  240.         return *this;
  241.     }
  242. }
  243.  
  244. /************************************************
  245.  *   --- sequence sequence::operator-  ---
  246.  *    Subtract two sequences
  247.  ************************************************/
  248. sequence sequence::operator- (sequence& arg)
  249. {
  250.     int i,lower=min(lower_limit,arg.lower_limit);
  251.     int upper=max(upper_limit,arg.upper_limit);
  252.     sequence temp(lower,upper);
  253.  
  254.     for (i=lower;i<=upper;i++) {
  255.         temp.stuff(i,fetch(i)-arg.fetch(i));
  256.     }
  257.     return temp;
  258. }
  259.  
  260. /************************************************
  261.  *  --- sequence& sequence::operator-=  ---
  262.  *    Subtract in place.  NOTE:  Call-by-reference
  263.  ************************************************/
  264. sequence& sequence::operator-= (sequence& arg)
  265. {
  266.     if ((lower_limit==arg.lower_limit)&&(upper_limit==arg.upper_limit)) {
  267.         int i,lim=upper_limit-lower_limit;
  268.         for (i=0;i<=lim;i++)seq[i]-=arg.seq[i];
  269.     } else {
  270.         sequence temp(*this);
  271.         *this=temp-arg;
  272.         return *this;
  273.     }
  274. }
  275.  
  276. /************************************************
  277.  *  --- sequence& sequence::operator=  ---
  278.  *    Copy operator.  NOTE, since we want separate copies of the
  279.  * signal, rather than just another POINTER to the SAME signal, we
  280.  * must allocate and copy all the data.
  281.  ************************************************/
  282. sequence& sequence::operator= (sequence& arg)
  283. {
  284.     int i;
  285.     /*we use old_seq, just in CASE somebody tries a=a*/
  286.     double huge *old_seq=seq;
  287.  
  288.     lower_limit=arg.lower_limit;
  289.     upper_limit=arg.upper_limit;
  290.     seq=new double[upper_limit+1-lower_limit];
  291.     if(seq==NULL){
  292.         cerr << "Cannot allocate sequence\n";
  293.         exit(1);
  294.     }
  295.     for(i=lower_limit;i<=upper_limit;i++){
  296.         stuff(i,arg.fetch(i));
  297.     }
  298.     delete old_seq;
  299.     return *this;
  300. }
  301.  
  302. /************************************************
  303.  *      --- int sequence::is_causal ---
  304.  *    Need this for certain stability tests, and IIR filtering
  305.  ************************************************/
  306. int sequence::is_causal(void)
  307. {
  308.     return(lower_limit>=0);
  309. }
  310.  
  311. sequence *fred, *mary, sam;
  312.  
  313. main()
  314. {
  315.     int u,l,i;
  316.  
  317.     fred= new sequence(0,5);
  318.     mary= new sequence(0,10);
  319.     for (i=0;i<=10;i++) {
  320.         fred->stuff(i,1.0);
  321.         mary->stuff(i,1.0);
  322.     }
  323.     cout << "\nconvolution:\n";
  324.     sam=(*fred)*(*mary);
  325.     sam.extent(&l,&u);
  326.     for (i=l;i<=u;i++) {
  327.         cout << "index=" << i << ", value=" <<sam.fetch(i) << "\n";
  328.     }
  329.     cout << "\naddition:\n";
  330.     sam=(*fred)+(*mary);
  331.     sam.extent(&l,&u);
  332.     for (i=l;i<=u;i++) {
  333.         cout << "index=" << i << ", value=" <<sam.fetch(i) << "\n";
  334.     }
  335.     cout << "\nautosubtration:\n";
  336.     sam-=(*fred)+(*mary);
  337.     sam.extent(&l,&u);
  338.     for (i=l;i<=u;i++) {
  339.         cout << "index=" << i << ", value=" <<sam.fetch(i) << "\n";
  340.     }
  341.     cout << "\nsubtraction:\n";
  342.     sam=(*fred)-(*mary);
  343.     sam.extent(&l,&u);
  344.     for (i=l;i<=u;i++) {
  345.         cout << "index=" << i << ", value=" <<sam.fetch(i) << "\n";
  346.     }
  347.     cout << "\nautoaddition:\n";
  348.     sam+=(*mary)-(*fred);
  349.     sam.extent(&l,&u);
  350.     for (i=l;i<=u;i++) {
  351.         cout << "index=" << i << ", value=" <<sam.fetch(i) << "\n";
  352.     }
  353.     delete fred;
  354.     delete mary;
  355.     fred= new sequence(0,200);
  356.     mary= new sequence(0,1);
  357.     fred->stuff(0,1);
  358.     mary->stuff(0,0.5);
  359.     mary->stuff(1,-0.495);
  360.     cout << "\ndivision:\n";
  361.     sam=(*fred)/(*mary);
  362.     sam.extent(&l,&u);
  363.     for (i=l;i<=u;i++) {
  364.         cout << "index=" << i << ", value=" <<sam.fetch(i) << "\n";
  365.     }
  366.     cout << "\nself division:\n";
  367.     (*fred)/=(*mary);
  368.     (*fred).extent(&l,&u);
  369.     for (i=l;i<=u;i++) {
  370.         cout << "index=" << i << ", value=" <<fred->fetch(i) << "\n";
  371.     }
  372.     cout << "\ndelta:\n";
  373.     (*fred)=sam-(*fred);
  374.     (*fred).extent(&l,&u);
  375.     for (i=l;i<=u;i++) {
  376.         cout << "index=" << i << ", value=" <<fred->fetch(i) << "\n";
  377.     }
  378.     cout << "making sequence\n";
  379.     delete fred;
  380.     fred=new sequence(sam);
  381.     cout << "\nbig multiplication\n";
  382.     sam=(*fred)*(*mary);
  383.     sam.extent(&l,&u);
  384.     for (i=l;i<=u;i++) {
  385.         cout << "index=" << i << ", value=" <<sam.fetch(i) << "\n";
  386.     }
  387.     sam=(*fred);
  388.     cout << "\n self multiplication\n";
  389.     sam*=(*mary);
  390.     sam.extent(&l,&u);
  391.     for (i=l;i<=u;i++) {
  392.         cout << "index=" << i << ", value=" <<sam.fetch(i) << "\n";
  393.     }
  394.     cout << "Ram left=" << farcoreleft() << "\n";
  395.     cout<<endl;
  396. }
  397.