diff --git a/Fortran/FUNC_sign.f90 b/Fortran/FUNC_sign.f90 new file mode 100644 index 0000000..19d7669 --- /dev/null +++ b/Fortran/FUNC_sign.f90 @@ -0,0 +1,19 @@ + + !符号函数 + real*8 function FUNC_sign(A) + implicit none + + real*8 :: A + if(isnan(A)) then + FUNC_sign = sqrt(-1.0) + elseif(A>0) then + FUNC_sign = 1.0 + elseif(A<0) then + FUNC_sign = -1.0 + elseif(A==0) then + FUNC_sign = 0.0 + end if + + end function + + \ No newline at end of file diff --git a/Fortran/SUB_WAVELET.f90 b/Fortran/SUB_WAVELET.f90 new file mode 100644 index 0000000..be871e7 --- /dev/null +++ b/Fortran/SUB_WAVELET.f90 @@ -0,0 +1,96 @@ + + !采用小波分解来自适应祛除随机误差 + 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 + + + + + \ No newline at end of file diff --git a/Fortran/SUB_appcoef.f90 b/Fortran/SUB_appcoef.f90 new file mode 100644 index 0000000..b356ccf --- /dev/null +++ b/Fortran/SUB_appcoef.f90 @@ -0,0 +1,115 @@ + + SUBROUTINE appcoef(nc,c,nl,l,n,IN3,na,a) + implicit none + + integer :: nc,nl,n,na + real*8 :: c(nc),a(1000) + integer :: l(nl) + character(len=3) :: IN3 + + integer :: rmax,nmax,nd,imax,lf,i + real*8 :: d(1000) + real*8,allocatable :: Lo_D(:),Hi_D(:),tempa(:) + + rmax = nl + nmax = rmax-2 + + na = 0 + a = 0 + + select case(IN3) + case('db2') + lf = 4 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.482962913144690,0.836516303737469,0.224143868041857,-0.129409522550921] + Hi_D = [-0.129409522550921,-0.224143868041857,0.836516303737469,-0.482962913144690] + + case('db3') + lf = 6 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.332670552950957,0.806891509313339,0.459877502119331,-0.135011020010391,-0.0854412738822415,0.0352262918821007] + Hi_D = [0.0352262918821007,0.0854412738822415,-0.135011020010391,-0.459877502119331,0.806891509313339,-0.332670552950957] + + case('db4') + lf = 8 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.230377813308855,0.714846570552542,0.630880767929590,-0.0279837694169839,-0.187034811718881,0.0308413818359870,& + 0.0328830116669829,-0.0105974017849973] + Hi_D = [-0.0105974017849973,-0.0328830116669829,0.0308413818359870,0.187034811718881,-0.0279837694169839,-0.630880767929590,& + 0.714846570552542,-0.230377813308855] + + case('db5') + lf = 10 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.160102397974125,0.603829269797473,0.724308528438574,0.138428145901103,-0.242294887066190,-0.0322448695850295,& + 0.0775714938400652,-0.00624149021301171,-0.0125807519990155,0.00333572528500155] + Hi_D = [0.00333572528500155,0.0125807519990155,-0.00624149021301171,-0.0775714938400652,-0.0322448695850295,0.242294887066190,& + 0.138428145901103,-0.724308528438574,0.603829269797473,-0.160102397974125] + + case('db6') + lf = 12 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.111540743350080,0.494623890398385,0.751133908021578,0.315250351709243,-0.226264693965169,-0.129766867567096,& + 0.0975016055870794,0.0275228655300163,-0.0315820393180312,0.000553842200993802,0.00477725751101065,-0.00107730108499558] + Hi_D = [-0.00107730108499558,-0.00477725751101065,0.000553842200993802,0.0315820393180312,0.0275228655300163,-0.0975016055870794,& + -0.129766867567096,0.226264693965169,0.315250351709243,-0.751133908021578,0.494623890398385,-0.111540743350080] + + case('db7') + lf = 14 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.0778520540850624,0.396539319482306,0.729132090846555,0.469782287405359,-0.143906003929106,-0.224036184994166,& + 0.0713092192670500,0.0806126091510659,-0.0380299369350346,-0.0165745416310156,0.0125509985560138,0.000429577973004703,& + -0.00180164070399983,0.000353713800001040] + Hi_D = [0.000353713800001040,0.00180164070399983,0.000429577973004703,-0.0125509985560138,-0.0165745416310156,0.0380299369350346,& + 0.0806126091510659,-0.0713092192670500,-0.224036184994166,0.143906003929106,0.469782287405359,-0.729132090846555,& + 0.396539319482306,-0.0778520540850624] + + case('db8') + lf = 16 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.0544158422430816,0.312871590914466,0.675630736298013,0.585354683654869,-0.0158291052560239,-0.284015542962428,& + 0.000472484573997973,0.128747426620186,-0.0173693010020221,-0.0440882539310647,0.0139810279170155,0.00874609404701566,& + -0.00487035299301066,-0.000391740372995977,0.000675449405998557,-0.000117476784002282] + Hi_D = [-0.000117476784002282,-0.000675449405998557,-0.000391740372995977,0.00487035299301066,0.00874609404701566,-0.0139810279170155,& + -0.0440882539310647,0.0173693010020221,0.128747426620186,-0.000472484573997973,-0.284015542962428,0.0158291052560239,& + 0.585354683654869,-0.675630736298013,0.312871590914466,-0.0544158422430816] + + case('fk4') + lf = 4 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.653927555569765,0.753272492839488,0.0531792287790598,-0.0461657148152177] + Hi_D = [-0.0461657148152177,-0.0531792287790598,0.753272492839488,-0.653927555569765] + + case('fk6') + lf = 6 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.427915032422310,0.812919643136907,0.356369511070187,-0.146438681272577,-0.0771777574069701,0.0406258144232379] + Hi_D = [0.0406258144232379,0.0771777574069701,-0.146438681272577,-0.356369511070187,0.812919643136907,-0.427915032422310] + + case('fk8') + lf = 8 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.349238111863800,0.782683620384065,0.475265135079471,-0.0996833284505732,-0.159978097434030,0.0431066681065162,& + 0.0425816316775818,-0.0190001788537359] + Hi_D = [-0.0190001788537359,-0.0425816316775818,0.0431066681065162,0.159978097434030,-0.0996833284505732,-0.475265135079471,& + 0.782683620384065,-0.349238111863800] + + end select + + na = l(1) + a(1:na) = c(1:na) + + imax = rmax+1 + do i=nmax,n+1,-1 + call detcoef(nc,c,nl,l,i,nd,d) + + allocate(tempa(l(imax-i))) + call idwt(na,a,nd,d,lf,Lo_D,Hi_D,l(imax-i),tempa) + na = l(imax-i) + a(1:na) = tempa + + deallocate(tempa) + end do + + deallocate(Lo_D,Hi_D) + end subroutine \ No newline at end of file diff --git a/Fortran/SUB_conv.f90 b/Fortran/SUB_conv.f90 new file mode 100644 index 0000000..ccd68dd --- /dev/null +++ b/Fortran/SUB_conv.f90 @@ -0,0 +1,18 @@ + + !向量卷积(全卷积) + SUBROUTINE conv(m,x,n,y,z) + implicit none + + integer :: m,n + real*8 :: x(m),y(n),z(m+n-1) + + integer :: i,j + + z = 0 + do i=1,m + do j=1,n + z(i+j-1) = z(i+j-1)+x(i)*y(j) + end do + end do + + end subroutine \ No newline at end of file diff --git a/Fortran/SUB_conv_valid.f90 b/Fortran/SUB_conv_valid.f90 new file mode 100644 index 0000000..14b6e68 --- /dev/null +++ b/Fortran/SUB_conv_valid.f90 @@ -0,0 +1,34 @@ + + !向量卷积(valid模式，去掉补零边缘) + !要求m>=n + !z的大小 = m-n+1 + SUBROUTINE conv_valid(m,x,n,y,z) + implicit none + + integer :: m,n + real*8 :: x(m),y(n),z(m-n+1) + + real*8 :: temp(m+n-1) + integer :: i,j + + !z = 0 + !do i=1,m-n+1 + ! do j=1,n + ! z(i) = z(i)+x(i+j-1)*y(j) + ! end do + !end do + + z = 0 + temp = 0 + do i=1,m + do j=1,n + temp(i+j-1) = temp(i+j-1)+x(i)*y(j) + end do + end do + + z = temp(n:m) + + end subroutine + + + \ No newline at end of file diff --git a/Fortran/SUB_cumsum.f90 b/Fortran/SUB_cumsum.f90 new file mode 100644 index 0000000..7455dba --- /dev/null +++ b/Fortran/SUB_cumsum.f90 @@ -0,0 +1,14 @@ + + !累积求和 + SUBROUTINE cumsum(n,A,B) + implicit none + + integer :: n,i + integer :: A(n),B(n) + + B(1) = A(1) + do i=2,n + B(i) = B(i-1)+A(i) + end do + + end subroutine \ No newline at end of file diff --git a/Fortran/SUB_detcoef.f90 b/Fortran/SUB_detcoef.f90 new file mode 100644 index 0000000..853e502 --- /dev/null +++ b/Fortran/SUB_detcoef.f90 @@ -0,0 +1,27 @@ + + SUBROUTINE detcoef(m,coefs,n,longs,levels,num,d) + implicit none + + integer :: m,n,levels,num + integer :: longs(n) + real*8 :: coefs(m),d(1000) + + integer :: nmax,tfirst(n),first(n-2),longs2(n-2),last(n-2),i + + d = 0 + + nmax = n-2 + call cumsum(n,longs,tfirst) + tfirst = tfirst+1 + + do i=n-2,1,-1 + first(n-i-1) = tfirst(i) + longs2(n-i-1) = longs(i+1) + end do + + last = first+longs2-1 + + num = last(levels)-first(levels)+1 + d(1:num) = coefs(first(levels):last(levels)) + + end subroutine \ No newline at end of file diff --git a/Fortran/SUB_idwt.f90 b/Fortran/SUB_idwt.f90 new file mode 100644 index 0000000..b4a0e3a --- /dev/null +++ b/Fortran/SUB_idwt.f90 @@ -0,0 +1,20 @@ + + SUBROUTINE idwt(na,a,nd,d,lf,Lo_D,Hi_D,L,x) + implicit none + + integer :: na,nd,lf,L + real*8 :: a(na),d(nd),Lo_D(lf),Hi_D(lf),x(L) + + real*8 :: temp1(L),temp2(L) + + x = 0 + temp1 = 0 + temp2 = 0 + + call upsconv1(na,a,lf,Lo_D,L,temp1) + call upsconv1(nd,d,lf,Hi_D,L,temp2) + + x = temp1+temp2 + + end subroutine + \ No newline at end of file diff --git a/Fortran/SUB_upsconv1.f90 b/Fortran/SUB_upsconv1.f90 new file mode 100644 index 0000000..aa3e1b5 --- /dev/null +++ b/Fortran/SUB_upsconv1.f90 @@ -0,0 +1,43 @@ + + SUBROUTINE upsconv1(nx,x,nf,f,s,y) + implicit none + + integer :: nx,nf,s + real*8 :: x(nx),f(nf),y(s) + + integer :: lx,lf,tempn,tempn2,i,first,last + real*8 :: d + real*8,allocatable :: tempy(:),tempy2(:) + + y = 0 + + lx = 2*nx + lf = nf + + tempn = lx-1 + allocate(tempy(tempn)) + tempy = 0 + do i=1,tempn,2 + tempy(i) = x(i/2+1) + end do + + tempn2 = tempn+nf-1 + allocate(tempy2(tempn2)) + call conv(tempn,tempy,nf,f,tempy2(1:tempn2)) + deallocate(tempy) + + if(s<0.or.s>=tempn2) then + y(1:tempn2) = tempy2(1:tempn2) + goto 10 + end if + + d = (tempn2-s)/2.0 + first = 1+int(floor(d)) + last = tempn2-int(ceiling(d)) + + !ny = last-first+1 + y(1:s) = tempy2(first:last) + + deallocate(tempy2) +10 end subroutine + \ No newline at end of file diff --git a/Fortran/SUB_wavedec.f90 b/Fortran/SUB_wavedec.f90 new file mode 100644 index 0000000..7978b55 --- /dev/null +++ b/Fortran/SUB_wavedec.f90 @@ -0,0 +1,141 @@ + + SUBROUTINE wavedec(s,x,n,IN3,nc,c,nl,l) + implicit none + + integer :: s,n,nc,nl + integer :: l(nl) + real*8 :: x(s),c(1000) + character(len=3) :: IN3 + real*8,allocatable :: Lo_D(:),Hi_D(:),y(:),a(:),d(:),z(:),temp(:) + integer :: lf,i,first,last,lenEXT,ny,nz,num,j,tempn + + c = 0 + l = 0 + + select case(IN3) + case('db2') + lf = 4 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [-0.129409522550921,0.224143868041857,0.836516303737469,0.482962913144690] + Hi_D = [-0.482962913144690,0.836516303737469,-0.224143868041857,-0.129409522550921] + + case('db3') + lf = 6 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.0352262918821007,-0.0854412738822415,-0.135011020010391,0.459877502119331,0.806891509313339,0.332670552950957] + Hi_D = [-0.332670552950957,0.806891509313339,-0.459877502119331,-0.135011020010391,0.0854412738822415,0.0352262918821007] + + case('db4') + lf = 8 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [-0.0105974017849973,0.0328830116669829,0.0308413818359870,-0.187034811718881,-0.0279837694169839,0.630880767929590,& + 0.714846570552542,0.230377813308855] + Hi_D = [-0.230377813308855,0.714846570552542,-0.630880767929590,-0.0279837694169839,0.187034811718881,0.0308413818359870,& + -0.0328830116669829,-0.0105974017849973] + + case('db5') + lf = 10 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.00333572528500155,-0.0125807519990155,-0.00624149021301171,0.0775714938400652,-0.0322448695850295,-0.242294887066190,& + 0.138428145901103,0.724308528438574,0.603829269797473,0.160102397974125] + Hi_D = [-0.160102397974125,0.603829269797473,-0.724308528438574,0.138428145901103,0.242294887066190,-0.0322448695850295,& + -0.0775714938400652,-0.00624149021301171,0.0125807519990155,0.00333572528500155] + + case('db6') + lf = 12 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [-0.00107730108499558,0.00477725751101065,0.000553842200993802,-0.0315820393180312,0.0275228655300163,0.0975016055870794,& + -0.129766867567096,-0.226264693965169,0.315250351709243,0.751133908021578,0.494623890398385,0.111540743350080] + Hi_D = [-0.111540743350080,0.494623890398385,-0.751133908021578,0.315250351709243,0.226264693965169,-0.129766867567096,& + -0.0975016055870794,0.0275228655300163,0.0315820393180312,0.000553842200993802,-0.00477725751101065,-0.00107730108499558] + + case('db7') + lf = 14 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.000353713800001040,-0.00180164070399983,0.000429577973004703,0.0125509985560138,-0.0165745416310156,-0.0380299369350346,& + 0.0806126091510659,0.0713092192670500,-0.224036184994166,-0.143906003929106,0.469782287405359,0.729132090846555,& + 0.396539319482306,0.0778520540850624] + Hi_D = [-0.0778520540850624,0.396539319482306,-0.729132090846555,0.469782287405359,0.143906003929106,-0.224036184994166,& + -0.0713092192670500,0.0806126091510659,0.0380299369350346,-0.0165745416310156,-0.0125509985560138,0.000429577973004703,& + 0.00180164070399983,0.000353713800001040] + + case('db8') + lf = 16 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [-0.000117476784002282,0.000675449405998557,-0.000391740372995977,-0.00487035299301066,0.00874609404701566,0.0139810279170155,& + -0.0440882539310647,-0.0173693010020221,0.128747426620186,0.000472484573997973,-0.284015542962428,-0.0158291052560239,& + 0.585354683654869,0.675630736298013,0.312871590914466,0.0544158422430816] + Hi_D = [-0.0544158422430816,0.312871590914466,-0.675630736298013,0.585354683654869,0.0158291052560239,-0.284015542962428,& + -0.000472484573997973,0.128747426620186,0.0173693010020221,-0.0440882539310647,-0.0139810279170155,0.00874609404701566,& + 0.00487035299301066,-0.000391740372995977,-0.000675449405998557,-0.000117476784002282] + + case('fk4') + lf = 4 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [-0.0461657148152177,0.0531792287790598,0.753272492839488,0.653927555569765] + Hi_D = [-0.653927555569765,0.753272492839488,-0.0531792287790598,-0.0461657148152177] + + case('fk6') + lf = 6 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [0.0406258144232379,-0.0771777574069701,-0.146438681272577,0.356369511070187,0.812919643136907,0.427915032422310] + Hi_D = [-0.427915032422310,0.812919643136907,-0.356369511070187,-0.146438681272577,0.0771777574069701,0.0406258144232379] + + case('fk8') + lf = 8 + allocate(Lo_D(lf),Hi_D(lf)) + Lo_D = [-0.0190001788537359,0.0425816316775818,0.0431066681065162,-0.159978097434030,-0.0996833284505732,0.475265135079471,& + 0.782683620384065,0.349238111863800] + Hi_D = [-0.349238111863800,0.782683620384065,-0.475265135079471,-0.0996833284505732,0.159978097434030,0.0431066681065162,& + -0.0425816316775818,-0.0190001788537359] + + end select + + l(n+2) = s + nc = 0 + tempn = s + allocate(temp(tempn)) + temp = x + do i=1,n + first = 2 + lenEXT = lf-1 + last = tempn+lf-1 + ny = tempn+2*lenEXT + nz = ny-lf+1 + allocate(y(ny),z(nz)) + call wextend(tempn,temp,lenEXT,y) + + call conv_valid(ny,y,lf,Lo_D,z) + num = (last-first)/2+1 + allocate(a(num),d(num)) + do j=1,num + a(j) = z(first+(j-1)*2) + end do + + call conv_valid(ny,y,lf,Hi_D,z) + do j=1,num + d(j) = z(first+(j-1)*2) + end do + + nc = num+nc + c(num+1:nc) = c(1:nc) + c(1:num) = d(1:num) + l(n+2-i) = num + + deallocate(temp) + tempn = num + allocate(temp(tempn)) + temp = a + + deallocate(a,y,z,d) + end do + + nc = num+nc + c(num+1:nc) = c(1:nc) + c(1:num) = temp + l(1) = num + + deallocate(temp,Lo_D,Hi_D) + end subroutine + + \ No newline at end of file diff --git a/Fortran/SUB_wextend.f90 b/Fortran/SUB_wextend.f90 new file mode 100644 index 0000000..31ec46a --- /dev/null +++ b/Fortran/SUB_wextend.f90 @@ -0,0 +1,28 @@ + + SUBROUTINE wextend(n,x,lf,y) + implicit none + + integer :: n,lf + real*8 :: x(n),y(n+2*lf) + + integer :: i + + y = 0 + + if(n>lf) then + do i=1,lf + y(i) = x(lf-i+1) + end do + do i=lf+1,lf+n + y(i) = x(i-lf) + end do + do i=lf+n+1,n+2*lf + y(i) = x(lf+2*n+1-i) + end do + else + goto 10 + end if + +10 end subroutine + + \ No newline at end of file diff --git a/Fortran/SUB_wthresh.f90 b/Fortran/SUB_wthresh.f90 new file mode 100644 index 0000000..c49b566 --- /dev/null +++ b/Fortran/SUB_wthresh.f90 @@ -0,0 +1,17 @@ + + SUBROUTINE wthresh(n,x,t,y) + implicit none + + integer :: n + real*8 :: x(n),t,y(n) + + real*8 :: temp(n),FUNC_sign + integer :: i + + temp = abs(x)-t + temp = (temp+abs(temp))/2.0 + do i=1,n + y(i) = FUNC_sign(x(i))*temp(i) + end do + + end subroutine \ No newline at end of file