#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\n",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