#include<stdio.h>
#include<math.h>
#include"mex.h"
#define Initial 1000
#define Pow2_Mem_Max 1024
#define Sur_Max 60
#define Pow2_Rate_Max 64
const double Gen_2[2]={369,491},Gen_3[3]={367,435,457};
void my_vit(double *In_Bit,double *Out_Sym,double *Gen,int Rate,int Mem_Len,int Sur_Len,int Fram_Len)
{
int i,j,k,m,sl,Xor_Self,Tmpr,flag;
static int Stat_Mac[Pow2_Mem_Max];/*[1<<Mem_Len]*/
int Tmp_Bit[Sur_Max][Pow2_Mem_Max];/*[Sur_Len*(1<<Mem_Len)]*/
double Min_Tmp;
double Tmp[Pow2_Rate_Max];/*[(1<<Rate)]*/
double Sum_Tmp1[Pow2_Mem_Max];/*[1<<Mem_Len]*/
static double Sum_Tmp[Pow2_Mem_Max];/*[(1<<Mem_Len)]*/
/*Initial part*/
for(i=0;i<(1<<Mem_Len);i++)
{Stat_Mac[i]=0;
for(j=0;j<Rate;j++)
{k=Xor_Self=0;
Tmpr = (i<<1)&(int)Gen[j];
while(Tmpr>>k)
Xor_Self ^= (Tmpr>>k++)&1;
Stat_Mac[i] += Xor_Self << j;
}
}
/*Processing*/
for(sl=0;sl<=Fram_Len/Sur_Len;sl++)
{
if((1+sl)*Sur_Len<=Fram_Len) m=Sur_Len;
else {m=Fram_Len%Sur_Len;if(m==0) break;}
for(i=0;i<m;i++)
{for(j=0;j<(1<<Rate);j++)
{Tmp[j]=0;
for(k= 0;k<Rate;k++)
Tmp[j] += ((((j>>k)&1)<<1)-1)*In_Bit[Rate*(i+sl*Sur_Len)+k];
}
for(k=0;k<(1<<Mem_Len);k++)
Sum_Tmp1[k]=Sum_Tmp[k];
for(k=0;k<((1<<Mem_Len)>>1);k++)
for(j=0;j<2;j++)
{ Tmp_Bit[i][(k<<1)+j] =(Sum_Tmp1[k]+Tmp[((1<<Rate)-1)*j-((j<<1)-1)*Stat_Mac[k]])<
(Sum_Tmp1[k+((1<<Mem_Len)>>1)]+Tmp[((1<<Rate)-1)*j-((j<<1)-1)*Stat_Mac[k+((1<<Mem_Len)>>1)]])?0:1;
Sum_Tmp[(k<<1)+j] = (Tmp_Bit[i][(k<<1)+j]==0)?(Sum_Tmp1[k]+Tmp[((1<<Rate)-1)*j-((j<<1)-1)*Stat_Mac[k]]):
(Sum_Tmp1[k+((1<<Mem_Len)>>1)]+Tmp[((1<<Rate)-1)*j-((j<<1)-1)*Stat_Mac[k+((1<<Mem_Len)>>1)]]);
}
}
Min_Tmp=Sum_Tmp[0];
for(k=1;k<(1<<Mem_Len);k++)
Min_Tmp=(Min_Tmp<=Sum_Tmp[k])?Min_Tmp:Sum_Tmp[k];
for(k=0;k<(1<<Mem_Len);k++)
if(Min_Tmp==Sum_Tmp[k])
{flag=k;
for(j=m-1;j>=0;j--)
{flag+=Tmp_Bit[j][flag]*(1<<Mem_Len);
if(sl*Sur_Len+j<Fram_Len-Mem_Len) Out_Sym[sl*Sur_Len+j] =(double)(flag&1);
flag >>= 1;
}
}
}
}
void mexFunction(int nlhs,mxArray *plhs[], int nrhs,const mxArray *prhs[])
{ double *Input,*Output,*Gen,Tmp,Tmp1[6];
double Mem_Len,Sur_Len;
int k,m,n,Gm,Gn,Chg;
Gen=Tmp1;
Input=mxGetPr(prhs[0]); /*(Length,Rate)*/
m=mxGetM(prhs[0]);
n=mxGetN(prhs[0]);
if(m>n) {Chg=n;n=m;m=Chg;}
Mem_Len=mxGetScalar(prhs[1]);
Sur_Len=mxGetScalar(prhs[2]);
Gm=mxGetM(prhs[3]);
Gn=mxGetN(prhs[3]);
if((Gm*Gn==1))
{Tmp=mxGetScalar(prhs[3]);
if(Tmp==2.0) for(k=0;k<2;k++) Tmp1[k]=Gen_2[k];
if(Tmp==3.0) for(k=0;k<3;k++) Tmp1[k]=Gen_3[k];
Gn=(int)Tmp;
}
else Gen=mxGetPr(prhs[3]);
if(Gm>Gn) Gn=Gm; /*Gn is Rate*/
if(nrhs<4)
mexErrMsgTxt("The number of input variable Must be at least four!\nYou shoud input in the following formation:\n convol(In_Symbol,Memory_Length,Gen_Polynomial)\nGood Luck! More detail,see convol.c");
plhs[0]= mxCreateDoubleMatrix(n/Gn-(int)Mem_Len,m,mxREAL);
Output=mxGetPr(plhs[0]);
for(k=0;k<m;k++)
my_vit(Input+k*n,Output+k*(n/Gn-(int)Mem_Len),Gen,Gn,(int)Mem_Len,(int)Sur_Len,n/Gn);
}