package signal.lib;

import java.util.Vector;

/**
   Same functionality as the hybrid class except that the parameters m, L, and EPS may be specified independently
   in each subinterval (local approach).  
*/
   


public final class hybridB {
  private static utility SF;
  private static chebyshevPolynomial C;
  private static gegenbauerPolynomial G;

//  In
// u      function passed to the method to be postprocessed
// d      location of edges (with d[0]=ltBoundaryPt and d[nd+1]=rtBoundaryPt), created by edgeDetect or manually input 
// mv     vector containing m in each subinterval which was created by hybrid.postProcess()
// Lv     vector containing L in each subinterval which was created by hybrid.postProcess()
// ev     vector containing EPS in each subinterval which was created by hybrid.postProcess()
// alpha  strength of the exponential filter exp(-alpha([k/N]^gamma))
// gamma  order of the exponential filter exp(-alpha([k/N]^gamma))

//   Out
// uG     the postprocessed function


   public static void postProcess(double u[], double d[], double uG[], Vector mv, Vector Lv,
                                        Vector ev, double alpha, double gamma) {
    int        N = u.length - 1;                                             
    double  ak[] = new double[N+1];   // Chebyshev coefficients
    double  aF[] = new double[N+1];   // Filtered Chebyshev Coefficients
    double   x[] = new double[N+1];
    int numDiscont = d.length;
    int sum=0;
    
     for ( int j=0; j<=N; j++) x[j] = C.xGLr(j,N);   // reverse order, i.e, 1.0...-1.0

// ---------------------------------------------------------------------------------   
// --- obtain Chebyshev coefficients for the function on the whole interval [-1,1] -
// ---------------------------------------------------------------------------------

      double dx=(1.0/N);
      C.chebyCoef(u,ak,N);
      
        double dxP=dx*Math.PI;
        double sum2=0;   

        for (int k=0; k<=N; k++ ) {
             sum2=0;
             for (int n=0; n<=N; n++ )  sum2 += C.CJ(n,N)*u[n]*Math.cos(k*n*dxP);
             aF[k] = 2.0*C.CJ(k,N)*dx*sum2;
          } // k
  
       for (int k=0; k<=N; k++)  aF[k] *= Math.exp(-alpha*Math.pow(k*dx,gamma));

// ---------------------------------------------------------------------------------
// ------  the reconstructed function  as  ug(i)= sum[ fh(L,m)*C(L,m,x) ] ----------
// ---------------------------------------------------------------------------------


      double L=0;  int m=0;
  

      int[] ma = new int[numDiscont];   double[] La = new double[numDiscont]; double[] ea = new double[numDiscont];
      Double Ld, ed; Integer md;
      for (int j=0; j<(numDiscont-1); j++) {
         md = (Integer)mv.elementAt(j);  Ld = (Double)Lv.elementAt(j);  ed = (Double)ev.elementAt(j); 
         ma[j] = md.intValue();          La[j] = Ld.doubleValue();      ea[j] = ed.doubleValue();
      } // j
      
      ed = (Double)ev.elementAt(numDiscont-1); 
      ea[numDiscont-1] = ed.doubleValue();
      
      double xi=0,K=0;
      
      int subinterval=0, oldSubInt=-1;
      for (int i=0; i<=N; i++) {
         uG[i]=0; xi=x[i];
         
          if (i==N) { // xi = -1.0
                       
                if (xi>=(d[0]+ea[0]) && xi<(d[1]-ea[1])) 
                  for (int k=0; k<=N; k++ )  uG[i] += aF[k]*Math.cos(i*k*dxP);
                
                   else {
            
                      subinterval=0; m = ma[0]; L = La[0];
                      if (subinterval!=oldSubInt) { oldSubInt=subinterval; K = G.km(m,L)/( G.km(m+1,L)*G.h(m,L) ); }

                      if (m==-1) uG[i] = u[i];                           // don't postprocess
                      else       uG[i] = gh(ak,m,L,N,d[0],d[1],xi,K);
            
                    } // else

         } else if (i==0) { // xi = 1.0
           
              if (xi>=(d[numDiscont-2]+ea[numDiscont-2]) && xi<(d[numDiscont-1]-ea[numDiscont-1])) 
                for (int k=0; k<=N; k++ )  uG[i] += aF[k]*Math.cos(i*k*dxP); 
                
                   else {
            
                      subinterval=numDiscont-2; m = ma[numDiscont-2]; L = La[numDiscont-2];
                      if (subinterval!=oldSubInt) { oldSubInt=subinterval; K = G.km(m,L)/( G.km(m+1,L)*G.h(m,L) ); }

                      if (m==-1) uG[i] = u[i];                           // don't postprocess
                      else       uG[i] = gh(ak,m,L,N,d[numDiscont-2],d[numDiscont-1],xi,K);
            
                    } // else
              
         } else {
         
             for (int j=0; j<(numDiscont-1); j++) {
               if (xi>=d[j] && xi<d[j+1]) {
         
                   if (xi>=(d[j]+ea[j]) && xi<(d[j+1]-ea[j+1]))  
                    for (int k=0; k<=N; k++ )  uG[i] += aF[k]*Math.cos(i*k*dxP);
                     //uG[i] = C.chebySeries(aF,xi);
                
                   else {
            
                      subinterval=j; m = ma[j]; L = La[j];
                      if (subinterval!=oldSubInt) { oldSubInt=subinterval; K = G.km(m,L)/( G.km(m+1,L)*G.h(m,L) ); }

                      if (m==-1) uG[i] = u[i];                           // don't postprocess
                      else       uG[i] = gh(ak,m,L,N,d[j],d[j+1],xi,K);
            
                    } // else
         
              break;    // get out of the j loop, the interval has been found
            } // if
         } // j
         } // else 

      } // i
    } // gegenReconstructHb
  
     

   private static double gh(double ak[], int m, double L, int N, double a, double b, double xi, double K) {
  
       int Nq = N + 2*(int)(L+1) + m;         
       double NI = 1.0/Nq;
       double sum=0, diff, gpL, g, gpm, gmmpL, gj, gpj, gp, xj, x;
       double eps = 0.5*(b-a), delta = 0.5*(b+a);
       
       x = G.tx(a,b,xi);
       
       for (int j=0; j<=Nq; j++) {
          xj = C.xGLr(j,Nq);
          diff = Math.abs(x-xj);
      
          if (diff<0.00000001) {    // if x == xj
          
                gpL  = NI*G.gegenPoly(m,L+1,xj);     g = Math.PI*K*G.gegenPoly(m,L,xj); 
                gpm  = NI*G.gegenPoly(m+1,L,xj);  gmmpL = Math.PI*K*G.gegenPoly(m-1,L+1,xj);
           
             sum += C.chebySeries(ak,eps*xj+delta)*Math.pow(1.0-xj*xj,L)*C.CJ(j,Nq)*2.0*L*(gpL*g-gpm*gmmpL);
      
          } else {
      
                gj  = NI*G.gegenPoly(m,L,xj);          gpj = NI*G.gegenPoly(m+1,L,xj); 
                 g  = Math.PI*K*G.gegenPoly(m,L,x);     gp = Math.PI*K*G.gegenPoly(m+1,L,x);
             
            sum += C.chebySeries(ak,eps*xj+delta)*Math.pow(1.0-xj*xj,L)*C.CJ(j,Nq)*( (gpj*g-gj*gp)/(xj-x) );
      
         } // if
      } // j
      
      return  sum;

  } // gh



} // hybridB