noise_process.c

#include <math.h>
#include "otdr_const.h"
#include "haar.h"
#include "otdr_bufferfilter.h"

WORD32 windowsize=200;
WORD32 backnum=900;
extern void filter(double *a, double *b, CurveData *data, int n);
void ReverseCurveData(CurveData *pCurveDataIn,CurveData *pCurveDataOut,WORD32 len)
{
WORD32 i;
i=len;
while(i--)
pCurveDataOut[len-(i+1)].y=pCurveDataIn[i].y;
}

void CurveMinusData(CurveData *pCurveDataIn,double data, WORD32 len)
{
WORD32 i=0;
for(i=0;i<len;i++)
pCurveDataIn[i].y=pCurveDataIn[i].y-data;
}

void FilterCurve(WORD32 plsWid, CurveData *pCurveData,const WORD32 len)
{
double a[2] = {0}, b[2] = {0};

NOISE_LPF_BUTTER(plsWid,a,b);
filter(a, b, pCurveData, len);
}

void CurveAddData(CurveData *pCurveDataIn,double data, WORD32 len)
{
WORD32 i=0;
for(i=0;i<len;i++)
pCurveDataIn[i].y=pCurveDataIn[i].y+data;
}

double getMean(CurveData *pCurveDataIn, WORD32 start,WORD32 n)
{
WORD32 i=0;
double ymean=0.0;
for(i=start;i<start+n;i++)
ymean+=pCurveDataIn[i].y/(double)n;
return ymean;
}

double getStd(CurveData *pCurveDataIn,WORD32 n)
{
WORD32 i=0;
double ymean=0.0,var=0.0,std=0.0;
for(;i<n;i++)
ymean+=pCurveDataIn[i].y/(double)n; /*getMean*/

for(i=0;i<n;i++)
var+=((pCurveDataIn[i].y-ymean)*(pCurveDataIn[i].y-ymean));

var/=(float)n;
std=sqrt(var);
return std;
}

WORD32 foundDataNo(CurveData *CurveDataIn,double ymean,double y2std,WORD32 found,WORD32 len)
{
int i=0;
double ymn=0.0;
for(i=0;i<len/windowsize-1;i++)
{
ymn=getMean(CurveDataIn,i*windowsize,windowsize);

if((ymn-ymean)>6*y2std)
{
if(i<(len/windowsize-1))
{
ymn=getMean(CurveDataIn,(i+1)*windowsize,windowsize);

if(ymn-ymean>3*y2std)
{
found=i;
break;
}
}
}
else if(ymn<ymean)
{
ymean=ymn;
}
}

return found;
}

void CaculateNoise(CurveData *pCurveDataIn1,CurveData *pCurveDataIn2,WORD32 found,WORD32 len,double *noise)
{
WORD32 j=0;
WORD32 k=found*windowsize-backnum ;
for(j=0;j<found*windowsize-backnum;j++)
{
noise[j]=pCurveDataIn1[len+backnum-found*windowsize+j].y - pCurveDataIn2[k--].y;
}
}

void ReconstructData(CurveData *pCurveDataIn,WORD32 found,WORD32 len,double *noise,double ymn,CurveData *pCurveDataOut)
{
WORD32 j = 0;
for(j = 0;j <=(len-found*windowsize+backnum-1);j++)
{
/* if(j>500 || j>found*windowsize-backnum)
break;*/
pCurveDataOut[j].y=pCurveDataIn[j].y-ymn;
}

for(j = (len-found*windowsize+backnum);j<len;j++)
{
pCurveDataOut[j].y=noise[j-(len-found*windowsize+backnum)];
}
for(j=0;j<len;j++)
{
pCurveDataIn[j].y=pCurveDataOut[j].y;
}
}

WORD32 NoiseProcess(Curve *pCurveOut,WORD32 sampinterval, WORD32 sampperiod, WORD32 len)
{
CurveData *data1=NULL;
CurveData *data2=NULL;
double ymean=0.0;
double y2std=0.0;
double ymn=0.0;
double *noise=NULL;
WORD32 found=0xff;
double endata = 0.0;

data1=(CurveData *)malloc(len * sizeof(CurveData));
if(NULL==data1)
{
return ERROR;
}
memset(data1,0,len);

data2=(CurveData *)malloc(len* sizeof(CurveData));
if(NULL==data2)
{
free(data1);
return ERROR;
}
memset(data2,0,len);

ReverseCurveData(pCurveOut->data,data1,len);
endata = getMean(pCurveOut,len-100,100);
CurveMinusData(data1,endata,len);

FilterCurve(pCurveOut->config.comm.plsWid,data1,len);

CurveAddData(data1,endata,len);

if(len<windowsize || len<500)
return ERROR;
ymean=getMean(data1,0,windowsize);
y2std=getStd(data1,500);

found=foundDataNo(data1,ymean,y2std,found,len);
printf("found:%d ",found);
if((found*windowsize-backnum-30)<0 || found == 0xff)
{
return ERROR;
}

ymn=getMean(data1,found*windowsize-backnum-30,30);

noise=(double *)malloc((found*windowsize-backnum)* sizeof(double));
if(NULL==noise)
{
free(data1);
free(data2);
return ERROR;
}
memset(noise,0,found*windowsize-backnum);

CaculateNoise(pCurveOut->data,data1,found,len,noise);

ReconstructData(pCurveOut->data,found,len,noise,ymn,data2);

free(data1);
free(data2);
free(noise);

return OK;

}
//没有专门的求任意底数对数的函数,不过可以用
//log(x)/log(y)表示log y x