C++编写S-G滤波

/*************************************************************************** * * Time: 2018-3-5 * Project: 滤波函数 * Purpose: 各种滤波函数的实现 * Author: wk * Copyright (c) 2018 * Describe: 实现各种滤波，以便后续的使用。目前仅实现了S-G滤波（该方法来源于文献Chen et al., 2004）， * 用SavitzkyGolay方法对EVI时序序列进行滤波去噪。 * ****************************************************************************/ #pragma once #include "Filter.h" Filter::Filter(){} Filter::~Filter(){} /** * @brief sg滤波 * @param y 平滑数组 * @param iSize 数组大小 * @return 返回值：表示计算过程中出现的各种错误信息 */ int Filter::SG_Filter(float *y,int iSize,float*Sgresult) { float *fQuality = new float[iSize]; memset(fQuality,0.0,sizeof(float)*iSize); findNoiceInEVI(y,iSize,0.2,fQuality); lineInterpol(y,iSize,fQuality); delete [] fQuality; //gsl_vector_free(gslQuality); //第三步：探测时序的长期趋势 //创建滤波系数 int iFilterSize = 2*5+1; float * fFilter = new float[iFilterSize]; memset(fFilter,0.0,sizeof(float)*iFilterSize); SG_FilterCreate(5,2,fFilter); //float * fFilter = new float[iFilterSize]; //for (int i= 0;i<iFilterSize;i++) // fFilter[i] = gsl_vector_get(sgFilter,i); //gsl_vector_free(sgFilter); //根据滤波进行卷积 float *ytr = new float[iSize]; memset(ytr,0.0,sizeof(float)*iSize); gsl_MyConv(y,iSize,fFilter,iFilterSize,ytr); delete [] fFilter; //cout<<"ytr0: "<<endl; //for (int i = 0;i<ytr->size;i++) //{ // cout<<gsl_vector_get(ytr,i)<<endl; //} //gsl_vector_free(sgFilter); //第四步：设定各点权重 //设定EVI序列各点的权重 //gsl_vector*gvWeight = gsl_CalWeight(y,ytr); float * fWeight = new float[iSize]; memset(fWeight,0.0,sizeof(float)*iSize); gsl_CalWeight(y,iSize,ytr,fWeight); //对初始结果的三次拟合结果进行判断，如果第一、二次的拟合指数最小，则滤波结束，否则进入滤波循环 ////////////////////////////////////////////////////////////////////////// float * ffittresulet = new float[3]; memset(ffittresulet,0.0,sizeof(float)*3); ////////////////////////////////////////////////////////////////////////// //第一次拟合 //第五步：产生一条新的EVI序列 float* newy1 = new float[iSize]; memset(newy1,0.0,sizeof(float)*iSize); gsl_GetMaxValueofVectors(y,iSize,ytr,newy1); //第六步：用SG滤波处理这条新的NDVI序列 //gsl_vector*sgFilter46 = SG_FilterCreate(4,6); int iFilter2Size = 2*4+1; float * fFilter2 = new float[iFilter2Size]; memset(fFilter2,0.0,sizeof(float)*iFilter2Size); SG_FilterCreate(4,6,fFilter2); //for(int i = 0 ; i<iFilter2Size;i++) // fFilter2[i] = gsl_vector_get(sgFilter46,i); //float *newy1 = new float[iSize]; gsl_MyConv(newy1,iSize,fFilter2,iFilter2Size,ytr); //gsl_vector_free(sgFilter46); //cout<<"ytr1: "<<endl; //for (int i = 0;i<ytr->size;i++) //{ // cout<<gsl_vector_get(ytr,i)<<endl; //} //第七步：计算拟合效果指数 ffittresulet[0]=gsl_CalTotal(ytr,iSize,y,fWeight); ////////////////////////////////////////////////////////////////////////// //第二次拟合 //第五步，产生一条新的EVI序列 float* newy2 = new float[iSize]; memset(newy2,0.0,sizeof(float)*iSize); gsl_GetMaxValueofVectors(y,iSize,ytr,newy2); //第六步：用SG滤波处理这条新的NDVI序列 gsl_MyConv(newy2,iSize,fFilter2,iFilter2Size,ytr); //cout<<"ytr2: "<<endl; //for (int i = 0;i<ytr->size;i++) //{ // cout<<gsl_vector_get(ytr,i)<<endl; //} //第七步：计算拟合效果指数 ffittresulet[1]= gsl_CalTotal(ytr,iSize,y,fWeight); ////////////////////////////////////////////////////////////////////////// //第三次拟合 float * newy3 = new float[iSize]; memset(newy3,0.0,sizeof(float)*iSize); gsl_GetMaxValueofVectors(y,iSize,ytr,newy3); gsl_MyConv(newy3,iSize,fFilter2,iFilter2Size,ytr); //cout<<"ytr3: "<<endl; //for (int i = 0;i<ytr->size;i++) //{ // cout<<gsl_vector_get(ytr,i)<<endl; //} ffittresulet[2]= gsl_CalTotal(ytr,iSize,y,fWeight); float minValue = findmin(ffittresulet,3); int iIndex = 0; for (int i = 0; i < 3; i++) { if (minValue == ffittresulet[i]) iIndex = i; } ////////////////////////////////////////////////////////////////////////// //用S-G方法进行循环滤波 float loop_times = 0; while(loop_times < 5 && iIndex!=1) { ffittresulet[0] = ffittresulet[1]; ffittresulet[1] = ffittresulet[2]; //gsl_vector_memcpy(newy1,newy2); //gsl_vector_memcpy(newy2,newy3); memoryCopy(newy2,newy1,iSize); memoryCopy(newy3,newy2,iSize); gsl_GetMaxValueofVectors(y,iSize,ytr,newy3); gsl_MyConv(newy3,iSize,fFilter2,iFilter2Size,ytr); //cout<<"ytr"<<4+loop_times<<": "<<endl; //for (int i = 0;i<ytr->size;i++) //{ // cout<<gsl_vector_get(ytr,i)<<endl; //} ffittresulet[2]= gsl_CalTotal(ytr,iSize,y,fWeight); //拟合效果指数f[1]是否最小？ float minValue = findmin(ffittresulet,3); for (int i = 0; i < 3; i++) { if (minValue == ffittresulet[i]) iIndex = i; } loop_times++; } for (int i =0;i<iSize;i++) Sgresult[i] = newy2[i]; delete [] ytr; delete [] newy1; delete [] newy2; delete [] newy3; delete [] fWeight; delete [] fFilter2; delete [] ffittresulet; //gsl_vector_free(ytr); //gsl_vector_free(sgFilter46); //gsl_vector_free(newy1); //gsl_vector_free(newy3); //gsl_vector_free(gvWeight); //cout<<"newy2"<<": "<<endl; //for (int i = 0;i<newy2->size;i++) //{ // cout<<gsl_vector_get(newy2,i)<<endl; //} return 0; } /** * 多项式拟合函数 * 根据 x y 拟合出一个N次多项式。返回多项式的系数。 */ int Filter::PolyFit(gsl_vector *x, gsl_vector *y, int N,float * a) { int i, j; int M = GSL_MIN(x->size, y->size); double chisq, xi; gsl_matrix *XX = gsl_matrix_alloc(M, N + 1); gsl_vector *c = gsl_vector_alloc(N + 1); gsl_matrix *cov = gsl_matrix_alloc(N + 1, N + 1); for(i = 0; i < M; i++) { gsl_matrix_set(XX, i, 0, 1.0); //获取gsl_vector*x的第i个值 xi = gsl_vector_get(x, i); for(j = 1; j <= N; j++) { gsl_matrix_set(XX, i, j, pow(xi, j)); } } gsl_multifit_linear_workspace *workspace = gsl_multifit_linear_alloc(M, N + 1); gsl_multifit_linear(XX, y, c, cov, &chisq, workspace); gsl_matrix_free(XX); gsl_matrix_free(cov); gsl_multifit_linear_free(workspace); for (int v =0;v<c->size;v++) a[v] = gsl_vector_get(c,v); gsl_vector_free(c); return 0; } // /** * 计算 Savitzky-Golay 滤波器系数 * SG 滤波器的阶数为 2M+1，多项式的最高次数为 N */ int Filter::SG_FilterCreate(int M, int N, float * sg) { int i; //创建2M+1的矩阵，并把元素全部初始化为0 gsl_vector *x = gsl_vector_alloc(2 * M + 1); gsl_vector *y = gsl_vector_alloc(2 * M + 1); //设置全是0 gsl_vector_set_zero(y); //设置y中第M个元素为1 gsl_vector_set(y, M, 1); for(i= -M; i <= M; i++) { gsl_vector_set(x, i + M, i); } float * c = new float[N+1]; memset(c,0.0,sizeof(float)*(N+1)); PolyFit(x, y, N,c); //释放 gsl_vector_free(x); gsl_vector_free(y); gsl_vector * b =gsl_vector_alloc(N+1); for (int j = 0;j<N+1;j++) gsl_vector_set(b,j,c[j]); delete [] c; gsl_vector *fir = gsl_vector_alloc(2 * M + 1); for(i = -M; i <= M; i++) { //gsl_poly_eval 多项式拟合，第一个为数组，第二个为多项式次数，第三个为自变量。c[N]*i^[N] gsl_vector_set(fir, i + M, gsl_poly_eval(b->data, N + 1, i)); } gsl_vector_free(b); for (int j=0;j<fir->size;j++) sg[j] = gsl_vector_get(fir,j); gsl_vector_free(fir); return 0; } //采用向两边无限扩展的方式处理边缘 // y = 1,2,3,4,5,6,7,8,9,10 // k = 1,2,3,4,5 //扩展y：1,1,///1,2,3,4,5,6,7,8,9,10,///10,10 //fArray,卷积数 //gvKernel.卷积核 //size : 卷积核大小 //fConv：结果存放 int Filter::gsl_MyConv(float*fArray,int iAsize,float*fKerne,int iKsize,float * fConv) { int addNum = (iKsize -1)/2; //creat a new expended y; //新向量gvNewy大小为: y->size+2*(k->size-1)/2 //gsl_vector *gvNewy = gsl_vector_alloc(size+2*addNum); float * fNewy = n