Convolução/Desconvolução de dados Fortran.
Caros colegas. Tenho aqui que executar um programa que após a leitura de dois ficheiros de dados, de 1024 linhas cada um, sendo um o correspondente aos dados experimentais e o outro o correspondente aos dados instrumentais, efectue a convolução ou a desconvolução dos mesmos. Encontrei uma série de numerical recepes sobre o assunto e construí um programa, cuja linha de código apresento abaixo, que lê correctamente os dados de entrada, mas no ficheiro de resultados criado, apresenta como valores para a desconvolução ou convolução, algo que se parece com tudo menos com o pretendido. Será que me poderão ajudar a encontrar o erro. Segue então de seguida o código do programa. Obrigado – Paulo Gomes
Caros colegas. Tenho aqui que executar um programa que após a leitura de dois ficheiros de dados, de 1024 linhas cada um, sendo um o correspondente aos dados experimentais e o outro o correspondente aos dados instrumentais, efectue a convolução ou a desconvolução dos mesmos. Encontrei uma série de numerical recepes sobre o assunto e construí um programa, cuja linha de código apresento abaixo, que lê correctamente os dados de entrada, mas no ficheiro de resultados criado, apresenta como valores para a desconvolução ou convolução, algo que se parece com tudo menos com o pretendido. Será que me poderão ajudar a encontrar o erro. Segue então de seguida o código do programa. Obrigado – Paulo Gomes
Código:
program glic
implicit none
CU USES convlv
integer cnt, i, w
character(10) mon,sct
integer m,n
real isign
real(8) :: scater,monomer
real(8) :: ans
real(8) :: data, respns
real(8) :: ress, maximo
dimension ress(2048)
dimension data(1024)
dimension respns(1024)
dimension ans(2048)
dimension scater(1024)
dimension monomer(1024)
common/data/mon,sct,scater,monomer,ress
cnt=1024
m=cnt
n=cnt
open(70,file='dados')
read(70,*) mon, sct
open (71, file=mon)
do i=1,cnt
read(71,*) monomer(i)
enddo
open (72, file=sct)
do i=1,cnt
read(72,*) scater(i)
enddo
call convlv(data,n,respns,m,isign,ans)
open(73, file='resultados')
do i=1,cnt
write(73,*) i, monomer(i), scater(i), ress(i)
enddo
end
SUBROUTINE convlv(data,n,respns,m,isign,ans)
INTEGER isign,m,n,NMAX
REAL data(n),respns(n)
COMPLEX ans(n)
character(10) mon,sct
real(8) :: scater, monomer
dimension scater(1024)
dimension monomer(1024)
real(8) :: ress
dimension ress(2048)
common/data/mon,sct,scater,monomer,ress
PARAMETER (NMAX=4096)
CU USES realft,twofft
INTEGER i,no2
COMPLEX fft(NMAX)
isign=-1
do i=1,1024
data(i)=monomer(i)
respns(i)=scater(i)
enddo
do 11 i=1,(m-1)/2
respns(n+1-i)=respns(m+1-i)
11 continue
do 12 i=(m+3)/2,n-(m-1)/2
respns(i)=0.0
12 continue
call twofft(data,respns,fft,ans,n)
no2=n/2
do 13 i=1,no2+1
if (isign.eq. 1) then
ans(i)=fft(i)*ans(i)/no2
ress(i)=ans(i)
else if (isign.eq.-1) then
if (abs(ans(i)).eq.0.0) pause
ans(i)=fft(i)/ans(i)/no2
ress(i)=ans(i)
else
pause 'no meaning for isign in convlv'
endif
13 continue
ans(1)=cmplx(real(ans(1)),real(ans(no2+1)))
ress(i)=ans(i)
call realft(ans,n,-1)
return
do i=1,1024
print(*,*)respns(i), data(i)
enddo
END
SUBROUTINE four1(data,nn,isign)
INTEGER isign,nn
REAL data(2*nn)
INTEGER i,istep,j,m,mmax,n
REAL tempi,tempr
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
n=2*nn
j=1
do 11 i=1,n,2
if(j.gt.i)then
tempr=data(j)
tempi=data(j+1)
data(j)=data(i)
data(j+1)=data(i+1)
data(i)=tempr
data(i+1)=tempi
endif
m=n/2
1 if ((m.ge.2).and.(j.gt.m)) then
j=j-m
m=m/2
goto 1
endif
j=j+m
11 continue
mmax=2
2 if (n.gt.mmax) then
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 13 m=1,mmax,2
do 12 i=m,n,istep
j=i+mmax
tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j)
data(j)=data(i)-tempr
data(j+1)=data(i+1)-tempi
data(i)=data(i)+tempr
data(i+1)=data(i+1)+tempi
12 continue
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
13 continue
mmax=istep
goto 2
endif
return
END
SUBROUTINE realft(data,n,isign)
INTEGER isign,n
REAL data(n)
CU USES four1
INTEGER i,i1,i2,i3,i4,n2p3
REAL c1,c2,h1i,h1r,h2i,h2r,wis,wrs
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
theta=3.141592653589793d0/dble(n/2)
c1=0.5
if (isign.eq.1) then
c2=-0.5
call four1(data,n/2,+1)
else
c2=0.5
theta=-theta
endif
wpr=-2.0d0*sin(0.5d0*theta)**2
wpi=sin(theta)
wr=1.0d0+wpr
wi=wpi
n2p3=n+3
do 11 i=2,n/4
i1=2*i-1
i2=i1+1
i3=n2p3-i2
i4=i3+1
wrs=sngl(wr)
wis=sngl(wi)
h1r=c1*(data(i1)+data(i3))
h1i=c1*(data(i2)-data(i4))
h2r=-c2*(data(i2)+data(i4))
h2i=c2*(data(i1)-data(i3))
data(i1)=h1r+wrs*h2r-wis*h2i
data(i2)=h1i+wrs*h2i+wis*h2r
data(i3)=h1r-wrs*h2r+wis*h2i
data(i4)=-h1i+wrs*h2i+wis*h2r
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
11 continue
if (isign.eq.1) then
h1r=data(1)
data(1)=h1r+data(2)
data(2)=h1r-data(2)
else
h1r=data(1)
data(1)=c1*(h1r+data(2))
data(2)=c1*(h1r-data(2))
call four1(data,n/2,-1)
endif
return
END
SUBROUTINE twofft(data1,data2,fft1,fft2,n)
INTEGER n
REAL data1(n),data2(n)
COMPLEX fft1(n),fft2(n)
CU USES four1
INTEGER j,n2
COMPLEX h1,h2,c1,c2
c1=cmplx(0.5,0.0)
c2=cmplx(0.0,-0.5)
do 11 j=1,n
fft1(j)=cmplx(data1(j),data2(j))
11 continue
call four1(fft1,n,1)
fft2(1)=cmplx(aimag(fft1(1)),0.0)
fft1(1)=cmplx(real(fft1(1)),0.0)
n2=n+2
do 12 j=2,n/2+1
h1=c1*(fft1(j)+conjg(fft1(n2-j)))
h2=c2*(fft1(j)-conjg(fft1(n2-j)))
fft1(j)=h1
fft1(n2-j)=conjg(h1)
fft2(j)=h2
fft2(n2-j)=conjg(h2)
12 continue
return
END
Última edição pelo moderador: