! ! FOURIER.f ! ! Fixed-Format Fortran Source File ! Generated by PGI Visual Fortran(R) ! 9/27/2013 9:05:33 AM ! program FOURIER implicit none integer i,k,Nint,Nmax,Nyquist,Ndati,iscratch,Nup,Ndown,j integer Nint_c integer Ndati_c integer itipo,iwin,ifilter,IOK parameter(Nmax=10001) integer iflg(Nmax) real Ak(Nmax/2),Bk(Nmax/2),Amp(Nmax/2),Fase(Nmax/2),filtro(Nmax) real t(Nmax),eta(Nmax),four(Nmax),scratch,alt(Nmax),scov real t_c(Nmax),eta_c(Nmax) real period,A0,omega,deltat,pigr,periodo_fisico,Energia real deltat_c,tbeg,tend character*30 filename write(*,*)'LEGGI DATI DA FILE (0) OPPURE GENERA DATI (1)?' read(*,*)itipo !************************************************************************** ! LEGGO IL FILE DATI E WINDOWING !************************************************************************** if(itipo.eq.0) then write(*,*)'DAMMI IL NOME DEL FILE CON ESTENSIONE' read(*,*)filename open(10,file=filename,status='old',access='sequential') i=1 99 read(10,*,end=66)iscratch i=i+1 go to 99 66 continue Ndati=i-1 write(*,*)'Ndati=',Ndati Nint=Ndati-1 if(Ndati.gt.Nmax) then write(*,*)'Ocio che il dimensionamento e sbagliato' stop end if rewind(10) do i=1,Ndati read(10,*)t(i),scratch,eta(i) end do close(10) open(10,file='serie_originale.dat',status='unknown', > access='sequential') do i=1,Ndati write(10,*)t(i),eta(i) end do close(10) Nint=Ndati-1 !-------------------------------------- write(*,*)'VUOI APPLICARE UNA SELEZIONE DEI DATI? (0=NO, 1=SI)' read(*,*)iwin if(iwin.eq.1) then write(*,*)'Dammi il tempo di inizio e di fine' read(*,*)tbeg,tend do i=1,Ndati iflg(i)=0 if(t(i).ge.tbeg.and.t(i).le.tend) iflg(i)=1 end do Ndati_c=0 do i=1,Ndati if(iflg(i).eq.1) then Ndati_c=Ndati_c+1 t_c(Ndati_c)=t(i) eta_c(Ndati_c)=eta(i) end if end do Ndati=Ndati_c do i=1,Ndati t(i)=t_c(i) eta(i)=eta_c(i) end do Nint=Ndati-1 open(10,file='serie_selezione.dat',status='unknown', > access='sequential') do i=1,Ndati write(10,*)t(i),eta(i) end do close(10) end if !************************************************************************** ! IN ALTERNATIVA GENERO INTERNAMENTE UNA SERIE TEMPORALE !************************************************************************** elseif(itipo.eq.1) then omega=1.0 pigr=acos(-1.) Periodo_fisico=2*pigr/omega period=Periodo_fisico*1.5 Nint_c=1000 deltat_c=period/float(Nint_c) Ndati_c=Nint_c+1 write(*,*)Ndati_c open(25,file='serie_generata.dat', > status='unknown',access='sequential') do i=1,Ndati_c t_c(i)=deltat_c*(i-1) eta_c(i)=5.0*cos(omega*t_c(i)) !SCRIVERE QUI LA FORMULA write(25,*)t_c(i),eta_c(i) end do close(25) Nint=Nint_c Ndati=Ndati_c do i=1,Ndati t(i)=t_c(i) eta(i)=eta_c(i) end do end if !************************************************************************** ! FINE LETTURA/GENERAZIONE/ DATI !************************************************************************** !************************************************************************** ! ANALISI IN FREQUENZA !************************************************************************** ! CALCOLO VARIABILI DI LAVORO period=t(Ndati)-t(1) omega=2*acos(-1.)/period write(*,*)'Periodo della REGISTRAZIONE= ',period ! CALCOLO A0 A0=0. do i=1,Nint A0=A0+(eta(i+1)+eta(i)) end do A0=A0/Nint !CALCOLO Ak Nyquist=Nint/2 do k=1,Nyquist Ak(k)=0.e0 Bk(k)=0.e0 do i=1,Nint Ak(k)=Ak(k)+(eta(i+1)*cos(k*omega*t(i+1))+ > eta(i)*cos(k*omega*t(i))) Bk(k)=Bk(k)+(eta(i+1)*sin(k*omega*t(i+1))+ > eta(i)*sin(k*omega*t(i))) end do Ak(k)=Ak(k)/Nint Bk(k)=Bk(k)/Nint Amp(k)=sqrt(Ak(k)**2+Bk(k)**2) Fase(k)=atan2(Bk(k),Ak(k)) end do !************************************************************************** ! FINE ANALISI IN FREQUENZA !************************************************************************** !************************************************************************** ! EVENTUALE FILTRATURA E RICOSTRUZIONE DEL SEGNALE FILTRATO !************************************************************************** 55 do k=1,Nyquist filtro(k)=1. end do write(*,*)'Indice di Nyquist=',Nyquist write(*,*) write(*,*)'APPLICAZIONE FILTRI (0=NO,1=SI)' read(*,*)ifilter if(ifilter.eq.1) then write(*,*)'Nup' read(*,*) Nup if(Nup.gt.Nyquist) Nup=Nyquist write(*,*)'Dammi Ndown' read(*,*) Ndown do k=1,Nyquist if(k.lt.Ndown) filtro(k)=0. if(k.gt.Nup) filtro(k)=0. end do open(10,file='filtro.dat', > status='unknown',access='sequential') do k=1,Nyquist write(10,*)k,k*omega,filtro(k) end do close(10) end if open(112,file='serie_Fourier.dat', > status='unknown',access='sequential') do i=1,Ndati four(i)=A0/2. do k=1,Nyquist four(i)=four(i)+filtro(k)*(Ak(k)*cos(k*omega*t(i))+ > Bk(k)*sin(k*omega*t(i))) end do write(112,'3(F15.5,1x)')t(i),four(i),eta(i) end do close(112) write(*,*)'SEI SODDISFATTO? 0=NO,1=SI' read(*,*)IOK if(IOK.eq.0) go to 55 101 continue !************************************************************************** ! FINE FILTRATURA !************************************************************************** !************************************************************************** ! CONSIDERAZIONI ENERGETICHE - IDENTITA' DI PARSEVAL !************************************************************************** Energia=0. do k=1,Nyquist Energia=Energia+Amp(k)**2 end do write(*,*)'Energia=',Energia !************************************************************************** ! STAMPA SU FILE LO SPETTRO DI AMPIEZZA E DI ENERGIA !************************************************************************** open(111,file='Fourier_AMP.dat', > status='unknown',access='sequential') write(111,'I4,1x,5(F15.5,1x)')0,omega*0,A0/2.,0.,A0/2.,0. do k=1,Nyquist write(111,'I4,1x,5(F15.5,1x)')k,omega*k,Ak(k),Bk(k), > Amp(k),Fase(k) end do close(111) stop open(111,file='Fourier_ENERGIA.dat', > status='unknown',access='sequential') write(111,'I4,1x,5(F15.5,1x)')0,omega*0,(A0/2.)**2, > (A0/2.)**2/omega do k=1,Nyquist write(111,'I4,1x,5(F15.5,1x)')k,omega*k,Amp(k)**2, > (Amp(k)**2)/omega end do close(111) !************************************************************************** ! !************************************************************************** stop end