#include <stdio.h>
#include <math.h>
#include "mex.h"

/****************************************************************/
/* CGTBANK                                                      */
/*         A complex-valued implementation of a filterbank      */
/*                                                              */
/*       By Aki Härmä April 13 1999     		        */
/****************************************************************/
void trans(double *ynr, double *rsignal, 
	   long int len,
	   double *AAr, double *AAi, int adim, double *BBr, 
	   double *BBi, int bdim,
	   double *rmem, double *imem, int N)
{
 
   int q, w, e, mlen, chan;
   long int o, aindex, bindex, outindex;
  double xr, xi, ffr, ffi, tmpr, tmpi;
  double Bbr,Bbi;
  double Ai[7],Ar[7],Br[7],Bi[7];

  aindex=bindex=outindex=0;
  if (adim>=bdim)mlen=adim; else mlen=bdim;
for (chan=0;chan<N;chan++)
  {
    for(q=0;q<adim;q++) {Ar[q]=AAr[aindex+N*q];Ai[q]=AAi[aindex+N*q];}
    for(q=0;q<bdim;q++) {Br[q]=BBr[bindex+N*q];Bi[q]=BBi[bindex+N*q];}
    aindex++;bindex++;
    for(q=0;q<mlen;q++) {rmem[q]=0;imem[q]=0;}
  Bbr=Br[0]/(Br[0]*Br[0]+Bi[0]*Bi[0]);
  Bbi=Bi[0]/(Br[0]*Br[0]+Bi[0]*Bi[0]);
   for(o=0;o<len;o++)
      {
       xr=rsignal[o];xi=0.0; ffr=0.0; ;
       /* update feedforward sum*/
       for(q=0; q<adim-1; q++) ffr+=Ar[q+1]*rmem[q]-Ai[q+1]*imem[q]; 
                                      
       /* update feedbackward sum*/
       for(q=0; q<bdim-1; q++) {xr-=Br[q+1]*rmem[q]-Bi[q+1]*imem[q]; 
                                xi-=Br[q+1]*imem[q]+ Bi[q+1]*rmem[q];}
      
       /* update output*/
       tmpr=xr;
       xr=Bbr*xr+Bbi*xi;xi=Bbr*xi-Bbi*tmpr;
       ynr[outindex++]=2*(Ar[0]*xr-Ai[0]*xi+ffr); 
       
       /* update inner states*/
       for(q=0;q<mlen-1;q++)
	 {
	   tmpr=rmem[q];tmpi=imem[q];
	   rmem[q]=xr;imem[q]=xi;
	   xr=tmpr;xi=tmpi;
	 }
      }
  }
return;
}

/************************  MEXFUNCTION *****************************/
void mexFunction(
	int		nlhs,
	mxArray	*plhs[],
	int		nrhs,
	const mxArray    *prhs[]
	)
{
	double	*yni, *ynr, *rsignal, *isignal;
	double *Ar,*Ai,*Br,*Bi;
	double *rmems,*imems;
	double *Rmems,*Imems;
        int q, n, m, adim, bdim, mlen, N;
        long int len,mm,nn;
/* GET POINTERS */
        Ar= mxGetPr(prhs[0]);
        Ai= mxGetPi(prhs[0]);
       	m = (int)mxGetM(prhs[0]);
	n = (int)mxGetN(prhs[0]);

	adim=n; N=m;
        Br= mxGetPr(prhs[1]);
        Bi= mxGetPi(prhs[1]);
      	m = (int)mxGetM(prhs[1]);
	n = (int)mxGetN(prhs[1]);
	bdim=n;
                    if(adim>=bdim)mlen=adim; else mlen=bdim;
	  rsignal = mxGetPr(prhs[2]);
	mm = (long int)mxGetM(prhs[2]);
	nn= (long int)mxGetN(prhs[2]);
	if(nn==1)len=mm; else len=nn;
       plhs[0] = mxCreateDoubleMatrix(len, N, mxREAL);
	ynr = mxGetPr(plhs[0]);
       plhs[1]= mxCreateDoubleMatrix(1,mlen+1,mxCOMPLEX);
	Rmems=mxGetPr(plhs[1]);
	Imems=mxGetPi(plhs[1]);

	/* Do the actual computations in a subroutine */
       if (nrhs==4){ 
	 rmems=mxGetPr(prhs[3]); imems= mxGetPi(prhs[3]);
	 for(q=0;q<mlen;q++) {Rmems[q]=rmems[q]; Imems[q]=imems[q];}
                          }
	trans(ynr,rsignal,len,Ar,Ai,adim,Br,Bi,bdim,Rmems,Imems,N);
	return;
}