96 lines
2.7 KiB
Fortran
96 lines
2.7 KiB
Fortran
|
|
|||
|
|
!<21><><EFBFBD><EFBFBD>С<EFBFBD><D0A1><EFBFBD>ֽ<EFBFBD><D6BD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><D3A6><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
SUBROUTINE SUB_WAVELET(n,&
|
|||
|
|
Data_Error,&
|
|||
|
|
Data_Re,&
|
|||
|
|
Error)&
|
|||
|
|
bind(C, name="SUB_WAVELET")
|
|||
|
|
!DEC$ ATTRIBUTES DLLEXPORT::SUB_WAVELET
|
|||
|
|
implicit none
|
|||
|
|
|
|||
|
|
!<21>ӿ<EFBFBD>
|
|||
|
|
integer :: n
|
|||
|
|
real*8 :: Data_Error(n),Data_Re(n),Error(n)
|
|||
|
|
|
|||
|
|
!<21>м<EFBFBD><D0BC><EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
character(len=3) :: Wavelets(10) !<21>˲<EFBFBD><CBB2><EFBFBD>ʽ
|
|||
|
|
integer :: Plies_Numbers(5) !<21>ֲ<EFBFBD><D6B2><EFBFBD>
|
|||
|
|
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)
|
|||
|
|
|
|||
|
|
!<21><>ʼ<EFBFBD><CABC>
|
|||
|
|
Data_Re = 0
|
|||
|
|
Error = 0
|
|||
|
|
|
|||
|
|
!<21><><EFBFBD><EFBFBD>
|
|||
|
|
!<21><><EFBFBD><EFBFBD><EFBFBD>źŷֽ<C5B7><D6BD><EFBFBD><EFBFBD>˲<EFBFBD><CBB2><EFBFBD>ʽ<EFBFBD>Լ<EFBFBD><D4BC>ֲ<EFBFBD><D6B2><EFBFBD>Ӱ<EFBFBD>죬<EFBFBD>˴<EFBFBD><CBB4><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ѱ<EFBFBD><D1B0>
|
|||
|
|
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
|
|||
|
|
!<21>źŷֽ<C5B7>
|
|||
|
|
nl = Plies_Numbers(j)+2
|
|||
|
|
allocate(l(nl))
|
|||
|
|
call wavedec(n,Data_Error,Plies_Numbers(j),Wavelets(i),nc,c,nl,l)
|
|||
|
|
|
|||
|
|
!<21><>ȡ<EFBFBD><C8A1><EFBFBD><EFBFBD>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><CFB5>
|
|||
|
|
call appcoef(nc,c(1:nc),nl,l,Plies_Numbers(j),Wavelets(i),nca,ca)
|
|||
|
|
|
|||
|
|
!<21><>ֵ<EFBFBD><D6B5>ȡ
|
|||
|
|
thr = (2.0*log(real(n)))**0.5
|
|||
|
|
|
|||
|
|
!<21><>ȡ<EFBFBD><C8A1><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ƶϸ<C6B5><CFB8>ϵ<EFBFBD><CFB5>
|
|||
|
|
ncd = 0
|
|||
|
|
do k=1,Plies_Numbers(j)
|
|||
|
|
call detcoef(nc,c,nl,l,k,ncd(k),cd(:,k))
|
|||
|
|
|
|||
|
|
!<21><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><D6B5><EFBFBD><EFBFBD>
|
|||
|
|
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
|
|||
|
|
|
|||
|
|
!<21><><EFBFBD><EFBFBD><EFBFBD>ع<EFBFBD>
|
|||
|
|
call appcoef(nc1,c1(1:nc1),nl,l,0,Wavelets(i),n,Data_Re)
|
|||
|
|
|
|||
|
|
!<21><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
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
|
|||
|
|
|
|||
|
|
!<21><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
Data_Re = Wave_Bset
|
|||
|
|
Error = Data_Error-Data_Re
|
|||
|
|
|
|||
|
|
end subroutine
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|