mathtools/windowdftcode.f90 at main · jjomiste/mathtools · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
program windowdft
  use mtmodules
  implicit none

  character(len=100):: infile
  integer:: ifile,xcol,ycol,outfile, wdftfile
  real(8), allocatable :: x(:), y(:), omega(:), fre(:), fim(:)
  real(8), allocatable :: ywindow(:)
  integer:: ij, jk, ik,tstep
  real(8) :: windl!! duration of the window
  real(8) :: peak !! peak of the Hann window
  real(8) :: xx !! variable
  character(len=50) :: outfilechar

  call getarg(1,infile) !! get the input file name

  !! set the columns to analyze
  write(*,*) 'Which is the column for the x?'
  read(*,*) xcol

  write(*,*) 'Which is the column for the y?'
  read(*,*) ycol

  write(*,*) 'What is the duration of the window?'
  read(*,*) windl

  write(*,*) 'Time steps in the to be stored'
  read(*,*) tstep

  !!avoiding errors in the duration of the window
  if (windl.le.1d-3) then
     write(*,*)
     write(*,*) 'THE DURATION OF THE PULSE MUST BE LARGER THAN ZERO'
     write(*,*) 'Breaking the code...'
     write(*,*)
     stop
  end if

  write(*,*) ''

  !!READING THE FILE
  call reading2list(infile,xcol,ycol,x,y)

  !!COMPUTING THE Windowed-DFT

  !! Open the new file
  write(outfilechar,'(A17,E9.3)') 'salida_wdft_file_',windl

  open(newunit=wdftfile, file=trim(outfilechar))

  write(wdftfile,*) '#Input file: ', trim(infile)
  write(wdftfile,*) ''
  write(wdftfile,*) '#Duration of the Hann Window: ', windl
  write(wdftfile,*) ''

  write(wdftfile,*) '# Note that the frequency units are the inverse of the input'
  write(wdftfile,'(A2,3X, A1,17X,A1,16X,A13,3X,A7,10X,A7,10X,A9)') ' #',  't','f', 'omega(2*pi*f)',  &
       'Re(FFT)',    'Im(FFT)',    '|WFFT|^2'
  write(wdftfile,*) '# '

  write(*,*) 'Computing the Windowed DFT'
  write(*,*)

  !! LOOP IN THE TIME POINTS
  Do ij=1,size(x),tstep

     peak=x(ij)
     !! Obtaining the time signal times the window
     call hannarray(x,y,peak,windl,ywindow)

     !! Compute the Fourier Transform of the tiem signal times the window
     call compute_wdft(x, ywindow,peak-windl*half,peak+windl*half, omega, fre, fim)

     do jk=1,size(omega)

        write(wdftfile,'(6E17.5)') x(ij), omega(jk)/(two*pi),omega(jk), fre(jk),fim(jk),fre(jk)**2.+fim(jk)**2.
     end do

     write(wdftfile,*) ''
  End Do


end program windowdft