//=============================================================
//     Numerical Analysis Packages in Java
//     Savitsky - Golay smoothing filter
//     Author : Dr. Turhan Coban
//=============================================================
import java.io.*;
import java.text.*;
import java.util.Locale;

public class savitsky_golay
{


public static double[] savgol(int np, int nl, int nr, int ld, int m)
	int NRANSI;
	int imj,ipj,j,k,kk,mm,*indx;
	float d,fac,sum,**a,*b;

	if (np < nl+nr+1 || nl < 0 || nr < 0 || ld > m || nl+nr < m)
	nrerror("bad args in savgol");
	int indx[]=new int[m+1];
	indx=ivector(1,m+1);
	a=matrix(1,m+1,1,m+1);
	b=vector(1,m+1);
	for (ipj=0;ipj<=(m << 1);ipj++) {
		sum=(ipj ? 0.0 : 1.0);
		for (k=1;k<=nr;k++) sum += pow((double)k,(double)ipj);
		for (k=1;k<=nl;k++) sum += pow((double)-k,(double)ipj);
		mm=FMIN(ipj,2*m-ipj);
		for (imj = -mm;imj<=mm;imj+=2) a[1+(ipj+imj)/2][1+(ipj-imj)/2]=sum;
	}
	ludcmp(a,m+1,indx,&d);
	for (j=1;j<=m+1;j++) b[j]=0.0;
	b[ld+1]=1.0;
	lubksb(a,m+1,indx,b);
	for (kk=1;kk<=np;kk++) c[kk]=0.0;
	for (k = -nl;k<=nr;k++) {
		sum=b[1];
		fac=1.0;
		for (mm=1;mm<=m;mm++) sum += b[mm+1]*(fac *= k);
		kk=((np-k) % np)+1;
		c[kk]=sum;
	}
	free_vector(b,1,m+1);
	free_matrix(a,1,m+1,1,m+1);
	free_ivector(indx,1,m+1);
}
	
//LU decomposition method
public static double[][] LU(double c[][],int indx[],int d[])
{
//returns LU decomposition of matrix c and index indx
double a[][];
int n=c.length;
a=new double[n][n];
double vv[];
vv=new double[n];
double sum,dum,big,temp;
int i,j,k;
int imax;
int nmax=100;
double tiny=1.0e-40;
imax=0;
for(i=1;i<=n;i++)
{
  for(j=1;j<=n;j++)
     a[i-1][j-1]=c[i-1][j-1];
}
d[0]=1;
for(i=1;i<=n;i++)
  {
  big=0.0;
  for(j=1;j<=n;j++)
    {
    if(Math.abs(a[i-1][j-1])>big) big=Math.abs(a[i-1][j-1]);
    }
  if(big==0) {System.out.println("singular matrix");return a;}
  vv[i-1]=1.0/big;
  }
for(j=1;j<=n;j++)
{
  for(i=1;i<j;i++)
    {
      sum=a[i-1][j-1];
      for(k=1;k<i;k++)
	{
	  sum-=a[i-1][k-1]*a[k-1][j-1];
	}
    a[i-1][j-1]=sum;
    }
  big=0;
  for(i=j;i<=n;i++)
    {
    sum=a[i-1][j-1];
    for(k=1;k<j;k++)
      {
      sum-=a[i-1][k-1]*a[k-1][j-1];
      }
    a[i-1][j-1]=sum;
    dum=vv[i-1]*Math.abs(sum);
    if(dum>=big)
      {
      imax=i;
      big=dum;
      }
    } //end of i=0
if(j != imax)
    {
    for(k=1;k<=n;k++)
      {
      dum=a[imax-1][k-1];
      a[imax-1][k-1]=a[j-1][k-1];
      a[j-1][k-1]=dum;
      }
    d[0]=-d[0];
    vv[imax-1]=vv[j-1];
    } //end of if
indx[j-1]=imax;
if(a[j-1][j-1]==0) a[j-1][j-1]=tiny;
if(j!=n)
  {
  dum=1.0/a[j-1][j-1];
  for(i=j+1;i<=n;i++)
    a[i-1][j-1]*=dum;
  }//endif
} //end for j=
return a;
}


public static double[] LUaxb(double a[][],double x[],int indx[])
{
//solves AX=B system of linear equation of LU decomposed matrix a
//(calculated by method LU)
int ii=0;
int i,j,ll=0;
double sum=0;
int n=a.length;
double b[];
b=new double[n];
for(i=1;i<=n;i++)
{
    b[i-1]=x[i-1];
}
for(i=1;i<=n;i++)
  {
  ll=indx[i-1];
  sum=b[ll-1];
  b[ll-1]=b[i-1];
  if(ii!=0)
     {
     for(j=ii;j<=(i-1);j++)
       {
       sum-=a[i-1][j-1]*b[j-1];
       }
     }
  else if(sum!=0) ii=i;
  b[i-1]=sum;
  }
for(i=n;i>=1;i--)
  {
  sum=b[i-1];
  if(i<n)
    {
    for(j=(i+1);j<=n;j++)
      {
      sum-=a[i-1][j-1]*b[j-1];
      }
    }
  b[i-1]=sum/a[i-1][i-1];
  }
return b;
}
}











