module fmft
implicit none
integer,parameter,private::NPMAX=2**12
integer,parameter,private::NPKMAX=20
integer,parameter,private::NPKMAX2=NPKMAX*2
double precision,private::dat0_mft(NPMAX),dat1_mft(NPMAX),dat_mft(NPMAX)
double precision,private::dat0_fmft1(NPMAX),dat1_fmft1(NPMAX),dat_fmft1(NPMAX)
double precision,private::dat0_fmft2(NPMAX),dat1_fmft2(NPMAX),dat_fmft2(NPMAX)
integer,private:: npts_mft,npts_fmft1,npts_fmft2 ! the real length of the input data
double precision,private:: dt_mft,dt_fmft1,dt_fmft2
double precision,private:: tt_mft,tt_fmft1,tt_fmft2
contains
subroutine mft(x,dt,pkfs,pkas,npk)
implicit none
integer::npk
integer::nh,ntry,nflag
integer:: i,n
double precision:: fr,fn ! spectrum reolution and the Nyquist frequency
double precision:: fmin,fmax,df,ftest ! boundary of the one-dimensional search and the stepsize
double precision:: ax,bx,cx,dx ! variables for one-dimensional search
double precision:: dt
double precision:: pkf,pka
double precision:: pkfs(:),pkas(:),x(:)
double precision:: pars(40)
call ini_mft(x,dt,nflag)
! print*,size(pkfs),npkmax
if(size(pkfs).gt.npkmax)then
nflag=0
endif
print*,"nflag=",nflag
call global_print()
fr=1d0/2d0/dt_mft/dble(npts_mft) ! resolution of the frequency
fn=1d0/2d0/dt_mft ! Nyquist frequency
npk=0
print*,"fr,fn",fr,fn
do while(npk.le.2)
npk=npk+1
call find_pkf_fft(dat_mft(1:npts_mft),dt,pkf)
fmax=pkf+fr
fmin=pkf-fr
df=(fmax-fmin)/1d2
if(fmax.ge.fn)fmax=fn-1d-10
if(fmin.le.0d0)fmin=1d-10
nh=1
do i=1,100
ftest=fmin+dble(i)*df
print*,ftest
print*,npts_mft,dt_mft
print*,ntry
print*,ftest,func_amp_hanning(ftest,pars,nh,ntry)
end do
enddo
end subroutine mft
subroutine ini_mft(x,dt,nflag)
implicit none
double precision:: x(:)
double precision:: dt,tt
integer::nflag,npts
nflag=1
dat0_mft(1:npmax)=0d0
dat_mft(1:npmax)=0d0
dat1_mft(1:npmax)=0d0
npts=size(x)
if(npts.gt.npmax)then
print*,"input error"
nflag=1
pause
return
endif
dt_mft=dt
dt_fmft1=dt
dt_fmft2=dt
npts_mft=npts
npts_fmft1=npts
npts_fmft2=npts
tt=dble(npts)*dt
tt_mft=tt
tt_fmft1=tt
tt_fmft2=tt
dat0_mft(1:npts)=x(1:npts)
dat_mft(1:npts)=dat0_mft(1:npts)
dat1_mft(1:npts)=dat0_mft(1:npts)
end subroutine ini_mft
subroutine find_pkf_fft(x,dt,pkf)
! find the peak frequency with the largest amplitude by employing the FFT(DFT) method
implicit none
integer:: i,i2
integer:: n,nhalf
double precision:: pkf,pka
double precision:: dt
double precision:: x(:)
double precision:: x0(size(x)),x1(size(x))
double precision:: freqs(size(x)/2),amps(size(x)/2)
n=size(x)
nhalf=n/2
x0(1:n)=x(1:n)
call drealft(x0,n,1) ! FFT to the frequency domain
x0(1:n)=x0(1:n)/dble(nhalf)
do i=1,nhalf-1
i2=i*2
freqs(i)=dble(i)/dble(n)/dt
amps(i)=dsqrt(x0(i2+1)**2+x0(i2+2)**2)
end do
freqs(nhalf)=1d0/dt/2d0
amps(nhalf)=dabs(x0(1)/2d0)
! do i=1,20
! print*,i,freqs(i),amps(i),x0(i),x(i)
! end do
pkf=freqs(2)
pka=amps(2)
do i=3,nhalf-1
if(amps(i).ge.pka)then
pka=amps(i)
pkf=freqs(i)
endif
end do
print*,"PKF,PKA",PKF,PKA
end subroutine find_pkf_fft
SUBROUTINE drealft(dat,n,isign)
INTEGER isign,n
DOUBLE PRECISION dat(n)
INTEGER i,i1,i2,i3,i4,n2p3
DOUBLE PRECISION c1,c2,h1i,h1r,h2i,h2r,wis,wrs
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
theta=3.141592653589793d0/dble(n/2)
c1=0.5d0
if (isign.eq.1) then
c2=-0.5d0
call dfour1(dat,n/2,+1)
else
c2=0.5d0
theta=-theta
endif
wpr=-2.0d0*sin(0.5d0*theta)**2
wpi=sin(theta)
wr=1.0d0+wpr
wi=wpi
n2p3=n+3
do i=2,n/4
i1=2*i-1
i2=i1+1
i3=n2p3-i2
i4=i3+1
wrs=wr
wis=wi
h1r=c1*(dat(i1)+dat(i3))
h1i=c1*(dat(i2)-dat(i4))
h2r=-c2*(dat(i2)+dat(i4))
h2i=c2*(dat(i1)-dat(i3))
dat(i1)=h1r+wrs*h2r-wis*h2i
dat(i2)=h1i+wrs*h2i+wis*h2r
dat(i3)=h1r-wrs*h2r+wis*h2i
dat(i4)=-h1i+wrs*h2i+wis*h2r
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
end do
if (isign.eq.1) then
h1r=dat(1)
dat(1)=h1r+dat(2)
dat(2)=h1r-dat(2)
else
h1r=dat(1)
dat(1)=c1*(h1r+dat(2))
dat(2)=c1*(h1r-dat(2))
call dfour1(dat,n/2,-1)
endif
return
END SUBROUTINE Drealft
SUBROUTINE dfour1(dat, nn, isign)
IMPLICIT NONE
INTEGER, INTENT(IN) :: nn, isign
DOUBLE PRECISION, DIMENSION(2*nn), INTENT(INOUT) :: dat
INTEGER :: i, j, m, mmax, n, istep
DOUBLE PRECISION :: tempr, tempi
DOUBLE PRECISION :: theta, wi, wpi, wpr, wr, wtemp
n = 2 * nn
j = 1
DO i = 1, n, 2
IF (j > i) THEN
tempr = dat(j)
tempi = dat(j+1)
dat(j) = dat(i)
dat(j+1) = dat(i+1)
dat(i) = tempr
dat(i+1) = tempi
END IF
m = n / 2
DO WHILE ((m >= 2) .AND. (j > m))
j = j - m
m = m / 2
END DO
j = j + m
END DO
mmax = 2
DO WHILE (n > mmax)
istep = 2 * mmax
theta = 6.28318530717959d0 / (isign * mmax)
wpr = -2.d0 * SIN(0.5d0 * theta)**2
wpi = SIN(theta)
wr = 1.d0
wi = 0.d0
DO m = 1, mmax, 2
DO i = m, n, istep
j = i + mmax
tempr = wr * dat(j) - wi * dat(j+1)
tempi = wr * dat(j+1) + wi * dat(j)
dat(j) = dat(i) - tempr
dat(j+1) = dat(i+1) - tempi
dat(i) = dat(i) + tempr
dat(i+1) = dat(i+1) + tempi
END DO
wtemp = wr
wr = wr * wpr - wi * wpi + wr
wi = wi * wpr + wtemp * wpi + wi
END DO
mmax = istep
END DO
END SUBROUTINE dfour1
subroutine global_print()
! subroutine used to print all the global variables
! defined inside the module
print*,"the initial length of the data is",NPMAX
print*,"the real length is",npts_mft
print*,"the sampling inteval is",dt_mft
print*,"the total time is",tt_mft
end subroutine global_print
function func_amp_hanning(fre,pars,nh,nflag)
implicit none
integer,intent(In):: nh
double precision,intent(in):: fre
integer:: i,nflag,ntry
double precision:: pi,dpi,rad
double precision:: t
double precision:: func_amp_hanning,amp,sf,cf
double precision:: pars(40)
print*,"inside the function func_amp_hanning"
pi=dacos(-1d0)
dpi=2d0*PI
rad=pi/180d0
nflag=1
! print*,"npts_mft in func_amp_hanning is",npts_mft
print*,"dat0 or dat,dat!!!"
! sf=0d0
! cf=0d0
! do i=1,1 !npts_mft
! t=dble(I)*dt_mft
! cf=cf+dat_mft(i)*dcos(dpi*fre*t)*(1d0-dcos(dpi*t/tt_mft))
! sf=sf+dat_mft(i)*dsin(dpi*fre*t)*(1d0-dcos(dpi*t/tt_mft))
! end do
! amp=dsqrt(cf**2+sf**2)
! pars(1)=cf
! pars(2)=sf
! func_amp_hanning=amp
return
end function func_amp_hanning
end module fmft
program main_fmft
use fmft
implicit none
integer,parameter::NPTS=2**10,NPTSH=NPTS/2,NPTSQ=NPTS/4
integer,parameter::NPKMAX=20,NPKMAX2=NPKMAX*2
integer:: i,ipk,npk
double precision::x0(NPTS),x(NPTS)
double precision::pi,dpi,rad
double precision:: dt,fre,t
double precision::pkf(NPKMAX),PKA(NPKMAX2)
pi=dacos(-1d0)
dpi=2d0*pi
rad=pi/180d0
dt=1d0
fre=1d0/2d0/dt/dble(nptsh)*10d0
do i=1,npts
t=dble(I)*dt
x0(i)=dcos(dpi*fre*t)
! print*,i,x0(i)
end do
print*,"fre is",fre
pause
call mft(x0,dt,pkf,pka,npk)
end program main_fmft
最新发布
本文介绍了一个使用汇编语言编写的简单程序，该程序从预定义的数据数组中查找并输出最小的偶数值。通过循环遍历数组并利用位操作判断每个元素是否为偶数，进而确定最小偶数。
扫一扫
2506
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
