Subject: スペクトル計算とハイパス＆ローパスフィルタ

sample.f90

! For compile
!  gfortran sample.f90 fftpack5.f90  -o sample -O3
!
!   need fftpack5.f90
!


program sample


 implicit none
 integer, parameter:: n=18000   ! data number
 integer, parameter:: ndiv=50  !   number of points in power&co-spectra after smoothing
 real :: x(n), y(n), dt=0.1
 real :: cospec(n), pow1(n), pow2(n), freq(n), azero1, azero2
 real :: cospec_smooth(ndiv), pow1_smooth(ndiv), pow2_smooth(ndiv), freq_smooth(ndiv)
 real :: x_high(n), cutoff, x_low(n)
 integer :: nh
 integer :: i, j 
 real, parameter :: pi = 3.14159265358979323846


! ****************************************************
! power spetra & co-spectrum

 call random_number(x)
 call random_number(y)

 

! do i=1, n
!  x(i) = 10*sin(2*pi/n*(3*i))
!  y(i) = 8*cos(2*pi/n*(5*i))
! end do

 call cospectrum( x, y, n, dt, cospec, pow1, pow2, freq, nh,  azero1, azero2 )

 open(unit=10, file="spectra_raw")
 do i=1, nh
  write(10,*) freq(i), pow1(i), pow2(i), cospec(i)
 end do
 close(10)

 call logsmooth( pow1, freq, nh, pow1_smooth, freq_smooth, ndiv)
 call logsmooth( pow2, freq, nh, pow2_smooth, freq_smooth, ndiv)
 call logsmooth( cospec, freq, nh, cospec_smooth, freq_smooth, ndiv)

 open(unit=10, file="spectra_smooth")
 do i=1, ndiv
  write(10,*) freq_smooth(i), pow1_smooth(i), pow2_smooth(i), cospec_smooth(i)
 end do
 close(10)


! ****************************************************
!  highpass filter

 cutoff = 0.5
 call random_number(x)
 call highpass(x, n, dt, cutoff, x_high )

 open(unit=10, file="highpass")
 do i=1, n
  write(10,*) x(i), x_high(i), x(i)-x_high(i)
 end do
close(10)


! ****************************************************
!  lowpass filter

 cutoff = 0.5
 ! call random_number(x)
 call lowpass(x, n, dt, cutoff, x_low )

 open(unit=10, file="lowpass")
 do i=1, n
  write(10,*) x(i), x_low(i), x(i)-x_low(i)
 end do
 close(10)

end

subroutine logsmooth( pow, freq, n, pow_smooth, freq_smooth, ndiv)
! smooth pow(n) into pow_smooth(ndiv)
!  same distance in log scale fot freq
!
implicit none
 real, parameter :: MISS = -9999
 integer :: n, ndiv
 real :: pow(n), freq(n)
 real :: pow_smooth(ndiv), freq_smooth(ndiv)
 real :: ldf, fmin, fmax
 integer :: i1, i, j, c
 real(kind=8) :: sf, sp

 !  determin fmin fmax
       ldf = ( log(freq(n)) - log(freq(1)) )/ndiv
       fmin = freq(1); i1=1

 do j=1, ndiv
  fmax = exp( ldf*j+log(freq(1)))
   sf = 0; sp = 0; c = 0; i = i1
!   write(0,*) "Average ", fmin, "-", fmax, "Hz"
   do while( freq(i) >= fmin .AND. freq(i) < fmax)
     sf = sf + freq(i); sp = sp + pow(i); c = c+1
     i = i + 1
   end do
   if( c > 0 ) then
     freq_smooth(j) = sf/dble(c); pow_smooth(j) = sp/dble(c)
   else
     freq_smooth(j) = MISS; pow_smooth(j) = MISS
   end if
  i1 = i; fmin = fmax
!   write(0,*) j,  freq_smooth(j),  pow_smooth(j)
 end do
end

 

 




subroutine lowpass(x, n, dt, cutoff, x_low )
implicit none
 integer:: n
 real :: x(n), x_low(n)
 real :: dt, cutoff

 real :: r(n)
 real ::  wsave(n+int(log(real(n)))+4), work(n)
 integer :: lenr, lensav, lenwork, i, j, nh, ier
 real :: freq

  lenr = n
  lensav = (n+int(log(real(n)))+4)
  lenwork = n


 call rfft1i( n, wsave, lensav, ier )
       if(ier /= 0) then
         write(0,*) 'ERR in rfft1i'
         stop
       end if

  r = x
  call rfft1f( n, 1, r, lenr, wsave, lensav, work, lenwork, ier )

! write(*,*) 'rfft1f'
         if( mod(n,2) == 0) then
             nh = n/2 -1
         else
             nh = (n-1)/2
         end if
!   write(*,*) 'n=', n, 'nh=', nh
 
         do i=1, nh
          freq = i/(n*dt)
!          write(*,*) "freq=", freq, "cutoff=", cutoff
          if( freq > cutoff ) then
              r(2*i-1+1) = 0   ! a
              r(2*i+1)   = 0   ! b
          end if
         end do
  
! write(*,*) 'go to rfft1b'
  call rfft1b( n, 1, r, lenr, wsave, lensav, work, lenwork, ier )
  x_low = r


end

subroutine highpass(x, n, dt, cutoff, x_high )
implicit none
 integer:: n
 real :: x(n), x_high(n)
 real :: dt, cutoff

 real :: r(n)
 real ::  wsave(n+int(log(real(n)))+4), work(n)
 integer :: lenr, lensav, lenwork, i, j, nh, ier
 real :: freq

  lenr = n
  lensav = (n+int(log(real(n)))+4)
  lenwork = n


 call rfft1i( n, wsave, lensav, ier )
       if(ier /= 0) then
         write(0,*) 'ERR in rfft1i'
         stop
       end if

  r = x
  call rfft1f( n, 1, r, lenr, wsave, lensav, work, lenwork, ier )

! write(*,*) 'rfft1f'
         if( mod(n,2) == 0) then
             nh = n/2 -1
         else
             nh = (n-1)/2
         end if
!   write(*,*) 'n=', n, 'nh=', nh
 
        r(1) = 0
         do i=1, nh
          freq = i/(n*dt)
!          write(*,*) "freq=", freq, "cutoff=", cutoff
          if( freq < cutoff ) then
              r(2*i-1+1) = 0   ! a
              r(2*i+1)   = 0   ! b
          end if
         end do
  
! write(*,*) 'go to rfft1b'
  call rfft1b( n, 1, r, lenr, wsave, lensav, work, lenwork, ier )
  x_high = r


end

subroutine  cospectrum( x, y, n, dt, cospec, pow1, pow2, freq, nh, azero1, azero2  )
implicit none

 integer:: n
 real :: x(n), y(n), cospec(n), pow1(n), pow2(n), freq(n), dt
 real :: azero1, azero2

 real :: r1(n), a1(n), b1(n)
 real :: r2(n), a2(n), b2(n)
 real :: wsave(n+int(log(real(n)))+4), work(n)
 integer :: lenr, lensav, lenwork, i, j, nh, ier

  lenr = n
  lensav = (n+int(log(real(n)))+4)
  lenwork = n


 call rfft1i( n, wsave, lensav, ier )
       if(ier /= 0) then
         write(0,*) 'ERR in rfft1i'
         stop
       end if

  r1 = x
 call rfft1f( n, 1, r1, lenr, wsave, lensav, work, lenwork, ier )
            if( ier /= 0 ) then
              write(0,*) 'ERROR in rfft1f', ier
              stop
             end if
 r2 = y;
 call rfft1f( n, 1, r2, lenr, wsave, lensav, work, lenwork, ier )
             if( ier /= 0 ) then
               write(0,*) 'ERROR in rfft1f', ier
               stop
              end if
	 azero1 = r1(1); azero2 = r2(1)
        if( mod(n,2) == 0) then
            nh = n/2 -1
        else
            nh = (n-1)/2
        end if

        do j=1,nh
         a1(j) = r1(2*j-1+1); a2(j) = r2(2*j-1+1)
         b1(j) = r1(2*j+1)  ; b2(j) = r2(2*j+1)
        end do

        cospec = (a1*a2+b1*b2)*(n*dt)
        pow1   = (a1*a1+b1*b1)*(n*dt)
        pow2   = (a2*a2+b2*b2)*(n*dt)
!   below just amplitude
!        cospec = (a1*a2+b1*b2)
!        pow1   = a1*a1+b1*b1
!        pow2   = a2*a2+b2*b2

        do i=1, nh
         freq(i) = i/(n*dt)
        end do
end