!采用小波分解来自适应祛除随机误差
    SUBROUTINE  SUB_WAVELET(n,&
							Data_Error,&
							Data_Re,&
							Error)&
	bind(C, name="SUB_WAVELET")
	!DEC$ ATTRIBUTES DLLEXPORT::SUB_WAVELET
    implicit none
    
    !接口
    integer :: n
    real*8 :: Data_Error(n),Data_Re(n),Error(n)
    
    !中间参数
    character(len=3) :: Wavelets(10)    !滤波形式
    integer :: Plies_Numbers(5)         !分层数
    real*8 :: c(1000),ca(1000),thr,cd(1000,6),ysoft(1000,6),c1(2000)
    integer :: nc,nl,nca,ncd(6),nc1,i,j,k,ii
    integer,allocatable :: l(:)
    real*8 :: R_SN,R_SN_Max
    character(len=3) :: Wavelets_best
    integer :: Plies_Bset(2)
    real*8 :: Wave_Bset(n)
    
    !初始化
    Data_Re = 0
    Error = 0
    
    !主体
    !由于信号分解受滤波形式以及分层的影响，此处进行穷举寻优
    Wavelets = ['db2','db3','db4','db5','db6',&
                'db7','db8','fk4','fk6','fk8']
    Plies_Numbers = [2,3,4,5,6]
    R_SN_Max = 0
    
    do i=1,10
        do j=1,5
            !信号分解
            nl = Plies_Numbers(j)+2
            allocate(l(nl))
            call wavedec(n,Data_Error,Plies_Numbers(j),Wavelets(i),nc,c,nl,l)
            
            !提取最后一层近似系数
            call appcoef(nc,c(1:nc),nl,l,Plies_Numbers(j),Wavelets(i),nca,ca)
            
            !阈值获取
            thr = (2.0*log(real(n)))**0.5
            
            !提取各层高频细节系数
            ncd = 0
            do k=1,Plies_Numbers(j)
                call detcoef(nc,c,nl,l,k,ncd(k),cd(:,k))
                
                !进行阈值处理
                call wthresh(ncd(k),cd(1:ncd(k),k),thr,ysoft(1:ncd(k),k))
            end do
            
            do k=1,Plies_Numbers(j)
                cd(:,k) = ysoft(:,k)
                nc1 = nca
                c1(1:nc1) = ca(1:nca)
                do ii=6,1,-1
                    if(ncd(ii)/=0) then
                        c1(nc1+1:nc1+ncd(ii)) = cd(1:ncd(ii),ii)
                        nc1 = nc1+ncd(ii)
                    end if
                end do
                
                !数据重构
                call appcoef(nc1,c1(1:nc1),nl,l,0,Wavelets(i),n,Data_Re)
                
                !计算信噪比
                R_SN = 10.0*log(sum(Data_Error**2)/sum((Data_Error-Data_Re)**2))
                if(R_SN > R_SN_Max) then
                    R_SN_Max = R_SN
                    Wavelets_best = Wavelets(i)
                    Plies_Bset(1) = Plies_Numbers(j)
                    Plies_Bset(2) = k
                    Wave_Bset = Data_Re
                end if
            end do
            deallocate(l)
        end do
    end do
    
    !计算误差      
    Data_Re = Wave_Bset
    Error = Data_Error-Data_Re
    
    end subroutine