! ftdemo - Discrete Fourier transform demonstration program
      program ftdemo
      integer*4 MAXN
      parameter( MAXN = 16384 )
      integer*4 N, i, j, k, method
      real*8 twoPiN, freq, phase, tau, pi, t(MAXN), y(MAXN)
      real*8 f(MAXN), powSpec(MAXN)
      complex*16 ii, yt(MAXN)

  !* Initialize the sine wave time series to be transformed.
      write(*,*) 'Enter the number of points: '
      read(*,*) N
      write(*,*) 'Enter frequency of the sine wave: '
      read(*,*) freq
      write(*,*) 'Enter phase of the sine wave: '
      read(*,*) phase
      tau = 1   ! Time increment
      pi = 3.141592654
      do i=1,N
        t(i) = (i-1)*tau                     ! t = [0, tau, 2*tau, ... ]
        y(i) = sin(2*pi*t(i)*freq + phase)   ! Sine wave time series
        f(i) = (i-1)/(N*tau)                 ! f = [0, 1/(N*tau), ... ]
      enddo

      !* Compute the transform using desired method: direct summation
      !  or fast Fourier transform (FFT) algorithm.
      write(*,*) 'Compute transform by, 1) Direct summation; 2) FFT: '
      read(*,*) method
      if( method .eq. 1 ) then            ! Direct summation
        twoPiN = -2*pi/N
        ii = (0., 1.)            ! ii = sqrt(-1)
        do k=0,N
          yt(k+1) = (0.0, 0.0)
          do j=0,(N-1)
            yt(k+1) = yt(k+1) +
     &                y(j+1)*cos(twoPiN*j*k) +
     &                ii*y(j+1)*sin(twoPiN*j*k)
          enddo
        enddo
      else                          ! Fast Fourier transform
        do i=1,N
          yt(i) = y(i)        ! Copy data for input to fft
        enddo
        call fft(yt,N)
      endif
      !* Compute the power spectrum
      do k=1,N
        powSpec(k) = abs(yt(k))**2
      enddo

      !* Print out the plotting variables:
      !    t, y, f, ytReal, ytImag, powspec
      open(11,file='t.txt',status='unknown')
      open(12,file='y.txt',status='unknown')
      open(13,file='f.txt',status='unknown')
      open(14,file='ytReal.txt',status='unknown')
      open(15,file='ytImag.txt',status='unknown')
      open(16,file='powSpec.txt',status='unknown')
      do i=1,N
        write(11,*) t(i)
        write(12,*) y(i)
        write(13,*) f(i)
        write(14,*) real(yt(i))
        write(15,*) imag(yt(i))
        write(16,*) powSpec(i)
      enddo
      stop
      end

!***** To plot in MATLAB; use the script below ********************
!load t.txt; load y.txt; load f.txt;
!load ytReal.txt; load ytImag.txt; load powSpec.txt;
!%* Graph the time series and its transform.
!figure(1); clf;  % Clear figure 1 window and bring forward
!plot(t,y);
!title('Original time series');
!ylabel('Amplitude');  xlabel('Time');
!figure(2); clf;  % Clear figure 2 window and bring forward
!plot(f,ytReal,'-',f,ytImag,'--');
!legend('Real','Imaginary');
!title('Fourier transform');
!ylabel('Transform');  xlabel('Frequency');
!%* Compute and graph the power spectrum of the time series.
!figure(3); clf;  % Clear figure 3 window and bring forward
!semilogy(f,powSpec,'-');
!title('Power spectrum (unnormalized)');
!ylabel('Power');  xlabel('Frequency');
!******************************************************************