misaki/SUB_WAVELET
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

修改算法

Browse Source
This commit is contained in:
misaki 2025-04-17 17:31:56 +08:00
parent 4c6342ff08
commit eae0c93c2f
12 changed files with 572 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
Fortran/FUNC_sign.f90 Normal file
View File

@ -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

96
Fortran/SUB_WAVELET.f90 Normal file
View File

@ -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

115
Fortran/SUB_appcoef.f90 Normal file
View File

@ -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

18
Fortran/SUB_conv.f90 Normal file
View File

@ -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

34
Fortran/SUB_conv_valid.f90 Normal file
View File

@ -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

14
Fortran/SUB_cumsum.f90 Normal file
View File

@ -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

27
Fortran/SUB_detcoef.f90 Normal file
View File

@ -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

20
Fortran/SUB_idwt.f90 Normal file
View File

@ -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

43
Fortran/SUB_upsconv1.f90 Normal file
View File

@ -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

141
Fortran/SUB_wavedec.f90 Normal file
View File

@ -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

28
Fortran/SUB_wextend.f90 Normal file
View File

@ -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

17
Fortran/SUB_wthresh.f90 Normal file
View File

@ -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