! ! 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_r,Nint_c integer Ndati_r,Ndati_c integer itipo parameter(Nmax=10001) real Ak(Nmax/2),Bk(Nmax/2),Amp(Nmax/2),Fase(Nmax/2) real t(Nmax),eta(Nmax),four(Nmax),scratch,alt(Nmax),scov real t_c(Nmax),eta_c(Nmax) real t_r(Nmax),eta_r(Nmax) real period,A0,omega,deltat,pigr,periodo_fisico,Energia real deltat_c,deltat_r,Periodo_ricampionamento write(*,*)'LEGGI DATI DA FILE (0) OPPURE GENERA DATI (1)?' read(*,*)itipo !************************************************************************** ! LEGGO IL FILE DATI !************************************************************************** if(itipo.eq.0) then ! open(10,file='Newyearwave2.txt',status='old',access='sequential') open(10,file='sonda.dat',status='old',access='sequential') ! open(10,file='gradino.dat',status='old',access='sequential') ! open(10,file='serie_generata.dat', ! > status='unknown',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(Nint+1.gt.Nmax) then write(*,*)'Ocio che il dimensionamento e sbagliato' stop end if rewind(10) do i=1,Nint+1 read(10,*)t(i),eta(i) end do close(10) !************************************************************************** ! 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*10.75 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)) 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 !************************************************************************** ! EVENTUALE RICAMPIONAMENTO DELLA SERIE SU UNA FINESTRA TEMPORALE OPPORTUNA !************************************************************************** go to 200 Nint_r=1000 Ndati_r=Nint_r+1 Periodo_ricampionamento = Periodo_fisico deltat_r=Periodo_ricampionamento/Nint_r do j=1,Nint_r+1 t_r(j)=t_c(1)+deltat_r*(j-1) end do do j=1,Nint_r+1 do i=1,Nint_c if(t_r(j).ge.t_c(i).and.t_r(j).le.t_c(i+1)) then eta_r(j)=eta_c(i)+ > (eta_c(i+1)-eta_c(i))/ > (t_c(i+1)-t_c(i))*(t_r(j)-t_c(i)) endif end do end do open(25,file='serie_ricampionata.dat', > status='unknown',access='sequential') do j=1,Ndati_r write(25,*)t_r(j),eta_r(j) end do close(25) Nint=Nint_r Ndati=Ndati_r do j=1,Ndati_r t(j)=t_r(j) eta(j)=eta_r(j) end do 200 continue !************************************************************************** ! FINE LETTURA/GENERAZIONE/RICAMPIONAMENTO DATI !************************************************************************** !************************************************************************** ! ANALISI IN FREQUENZA !************************************************************************** ! CALCOLO VARIABILI DI LAVORO period=t(Nint+1)-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 !************************************************************************** go to 100 55 write(*,*)'Dammi Nup' read(*,*) Nup if(Nup.gt.Nyquist) Nup=Nyquist write(*,*)'Dammi Ndown' read(*,*) Ndown open(112,file='Serie_Fourier.out', > status='unknown',access='sequential') do i=1,Nint+1 four(i)=A0/2. do k=Ndown, Nup !k=1,Nyquist ! write(*,*)k four(i)=four(i)+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) ! go to 55 100 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.out', > 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) open(111,file='Fourier_ENERGIA.out', > 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