1. Este site usa cookies. Ao continuar a usar este site está a concordar com o nosso uso de cookies. Saber Mais.

Convolução/Desconvolução de dados Fortran.

Discussão em 'Programação' iniciada por PGomes76, 29 de Janeiro de 2008. (Respostas: 0; Visualizações: 1495)

  1. 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


    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: 29 de Janeiro de 2008

Partilhar esta Página