package signal.lib; import java.util.Vector; /** The hyprid class contains one public method postProcess() which applies the Gegenbauer Reconstruction Procedure (GRP) in regions (d-EPS,d+EPS) and an exponential filter exp(-alpha([k/N]^gamma)) at points which fall outside these intervals. The reconstruction parameters are specified by a global approach, as in gegenbauerReconstruction. After this method is run on the data once, the hybridB.postProcess method may be run and the reconstruction parameters may be set separately in each interval (local approach). @see "A Hybrid Approach to Spectral Reconstruction of Piecewise Smooth Functions. Preprint 2000. Anne Gelb." */ public final class hybrid { 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 // LK Gegenbauer Exponent parameter in L = 0.5*LK*(d[j+1]-d[j])*N // mK parameter in order of the Gegenbauer Polynomial m = 0.5*MK*(d[j+1]-d[j])*N // alpha strength of the exponential filter exp(-alpha([k/N]^gamma)) // gamma order of the exponential filter exp(-alpha([k/N]^gamma)) // EPS the GRP is applied in intervals (d-EPS,d+EPS); in hybridB EPS may be specified differently in each subinterval // Out // uG the postprocessed data // mv vector of m in subintervals (mv can then be passed to hybridB as an initial choice of m in each subinterval) // Lv vector of L in subintervals (Lv vector can then be passed to hybridB as an initial choice of L in each subinterval) public static void postProcess(double u[], double d[], double LK, double MK, double uG[], Vector mv, Vector Lv, double alpha, double gamma, double EPS) { 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] - // --------------------------------------------------------------------------------- C.chebyCoef(u,ak,N); double dx=(1.0/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]; for (int j=0; j<(numDiscont-1); j++) { m = (int)Math.ceil(0.5*MK*(d[j+1]-d[j])*N); L = 0.5*LK*(d[j+1]-d[j])*N; mv.addElement(new Integer(m)); Lv.addElement(new Double(L)); ma[j]=m; La[j]=L; } // j 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]+EPS) && xi<(d[1]-EPS)) for (int k=0; k<=N; k++ ) uG[i] += aF[k]*Math.cos(i*k*dxP); // uG[i] = C.chebySeries(aF,xi); 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) ); } uG[i] = gh(ak,m,L,N,d[0],d[1],xi,K); } } else if (i==0) { // xi = 1.0 if (xi>=(d[numDiscont-2]+EPS) && xi<(d[numDiscont-1]-EPS)) for (int k=0; k<=N; k++ ) uG[i] += aF[k]*Math.cos(i*k*dxP); // uG[i] = C.chebySeries(aF,xi); 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) ); } uG[i] = gh(ak,m,L,N,d[numDiscont-2],d[numDiscont-1],xi,K); } } else { for (int j=0; j<(numDiscont-1); j++) { if (xi>=d[j] && xi=(d[j]+EPS) && xi<(d[j+1]-EPS)) 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) ); } uG[i] = gh(ak,m,L,N,d[j],d[j+1],xi,K); } break; // get out of the j loop, the interval has been found } // if } // j } // else } // i } // gegenReconstructH 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, xGL, 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 } // hybrid