96 lines
2.7 KiB
Fortran
96 lines
2.7 KiB
Fortran
|
|
!采用小波分解来自适应祛除随机误差
|
|
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
|
|
|
|
|
|
|
|
|
|
|