home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / dcosine.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1989-08-18  |  4.0 KB  |  168 lines
  1. /*
  2.  *    Definitions for the double precision cosine server
  3.  *
  4.  *    Copyright (C) 1988, 1989.
  5.  *
  6.  *    Dr. Thomas Keffer
  7.  *    Rogue Wave Associates
  8.  *    P.O. Box 85341
  9.  *    Seattle WA 98145-1341
  10.  *
  11.  *    Permission to use, copy, modify, and distribute this
  12.  *    software and its documentation for any purpose and
  13.  *    without fee is hereby granted, provided that the
  14.  *    above copyright notice appear in all copies and that
  15.  *    both that copyright notice and this permission notice
  16.  *    appear in supporting documentation.
  17.  *    
  18.  *    This software is provided "as is" without any
  19.  *    expressed or implied warranty.
  20.  *
  21.  *
  22.  *    @(#)dcosine.cc    2.1    8/18/89
  23.  */
  24.  
  25.  
  26.  
  27. #include "rw/DCosineFFT.h"
  28. #include "rw/mathpack.h"
  29.  
  30. DoubleCosineServer::DoubleCosineServer()
  31. {
  32.   server_N=0;
  33. }
  34.  
  35. DoubleCosineServer::DoubleCosineServer(unsigned N)
  36. {
  37.   setOrder(N);
  38. }
  39.  
  40. DoubleCosineServer::DoubleCosineServer(const DoubleCosineServer& s)
  41.      : (s), sins(s.sins)
  42. {
  43.   server_N = s.server_N;
  44. }
  45.  
  46. void
  47. DoubleCosineServer::operator=(const DoubleCosineServer& s)
  48. {
  49.   DoubleFFTServer::operator=(s);
  50.   server_N = s.server_N;
  51.   sins.reference(s.sins);
  52. }
  53.  
  54. void
  55. DoubleCosineServer::setOrder(unsigned N)
  56. {
  57.   server_N = N;
  58.   DoubleFFTServer::setOrder(N/2+1);
  59.   sins.reference(2.0*sin(DoubleVec(N-1,M_PI/N, M_PI/N)));
  60. }
  61.  
  62. /* 
  63. The fourier transform of a real even sequence is a cosine transform.
  64. The result is also a real even sequence.  This is prodedure #6 from
  65. Cooley, et al. (1970) J. Sound Vib. (12) 315--337.
  66. */
  67.  
  68. DoubleVec
  69. DoubleCosineServer::cosine(const DoubleVec& v)
  70. {
  71.   unsigned Nstore    = v.length();
  72.   unsigned N        = Nstore-1;
  73.   unsigned Nm1        = N-1;
  74.   unsigned Nhalf    = N/2;
  75.   unsigned Nhalfp1    = Nhalf+1;
  76.   unsigned Nhalfm1    = Nhalf-1;
  77.  
  78.   // Various things that could go wrong:
  79.   checkOdd(Nstore);
  80.   if( N != server_N )setOrder(N);
  81.  
  82.   /* The results will appear here, and will consist of length() points.*/
  83.   DoubleVec results(Nstore);
  84.  
  85.   /* X will be a complex conjugate even vector */
  86.   DComplexVec X(Nhalfp1);
  87.  
  88.   X.slice(1,Nhalfm1,1) =
  89.     DComplexVec(v.slice(2,Nhalfm1,2),    // Real part
  90.         v.slice(1,Nhalfm1,2) - v.slice(3,Nhalfm1,2)); // Imag part
  91.  
  92.   X(0) = DComplex(v(0));
  93.   X(Nhalf) = DComplex(v(N));
  94.  
  95.   /* Take the IDFT. The transform of a complex conjugate even vector
  96.     N/2+1 points long is a real vector, N points long. */
  97.   DoubleVec A = DoubleFFTServer::ifourier(X);
  98.   
  99.   DoubleVec temp1 = A.slice(1,Nm1,1);
  100.   DoubleVec temp2 = A.slice(Nm1,Nm1,-1);
  101.  
  102.   results.slice(1,Nm1,1) = 0.5*( (temp1 + temp2) - (temp1 - temp2)/sins);
  103.  
  104.   /* Special procedure required for the first and last points. */
  105.   // Sum all of the odd points.  Count all but v(N) twice.
  106.   double A2 = 2.0*sum(v.slice(1,Nhalf,2));
  107.  
  108.   if( N%2) A2 += v(N);        // If N is odd, include v(N)
  109.  
  110.   results(0) = A(0) + A2;
  111.   results(N) = A(0) - A2;
  112.   
  113.   return results;
  114. }
  115.  
  116. /* 
  117. The fourier transform of a real odd sequence is a sine transform.
  118. The result is also a real odd sequence.  This is prodedure #7 from
  119. Cooley, et al. (1970) J. Sound Vib. (12) 315--337.
  120. */
  121.  
  122. DoubleVec
  123. DoubleCosineServer::sine(const DoubleVec& v)
  124. {
  125.   unsigned Nstore    = v.length();
  126.   unsigned N        = Nstore+1;
  127.   unsigned Nm1        = N-1;
  128.   unsigned Nhalf    = N/2;
  129.   unsigned Nhalfp1    = Nhalf+1;
  130.   unsigned Nhalfm1    = Nhalf-1;
  131.  
  132.   // Various things that could go wrong:
  133.   checkOdd(Nstore);
  134.   if(N != server_N)setOrder(N);
  135.  
  136.   /* X will be a complex conjugate even vector */
  137.   DComplexVec X(Nhalfp1);
  138.  
  139.   X.slice(1,Nhalfm1,1) =
  140.     DComplexVec(v.slice(0,Nhalfm1,2) - v.slice(2,Nhalfm1,2),    // real part
  141.         -v.slice(1,Nhalfm1,2));                // imag part
  142.  
  143.   X(0)     = DComplex(-2.0*v(0));
  144.   X(Nhalf) = DComplex( 2.0*v(N-2));
  145.  
  146.   /* Take the IDFT. The transform of a complex conjugate even vector
  147.     N/2+1 points long is a real vector, N points long. */
  148.   DoubleVec A = DoubleFFTServer::ifourier(X);
  149.   
  150.   DoubleVec temp1 = A.slice(1,Nm1,1);
  151.   DoubleVec temp2 = A.slice(Nm1,Nm1,-1);
  152.  
  153.   DoubleVec results = 0.5*( (temp1 - temp2) - (temp1 + temp2)/sins);
  154.  
  155.   return results;
  156. }
  157.  
  158. void
  159. DoubleCosineServer::checkOdd(int l)
  160. {
  161.   if( !(l%2) ){
  162.     char msg[80];
  163.     sprintf(msg, "Length (%d) of a real series must be odd.", l);
  164.     RWnote("DoubleCosineServer::[i]fourier()", msg);
  165.     RWerror(DEFAULT);
  166.   }
  167. }
  168.