プログラムとサンプルファイル
Lanczos.tar.gz - 2011/12/23 20:20、Atsuyoshi Manda (バージョン 1) 削除
20KB ダウンロード
結果の例
lp_sf_eddy_data.ps - 2011/12/23 20:20、Atsuyoshi Manda (バージョン 1) 削除
29KB 表示 ダウンロード
使用例
[Sat Dec 24 10:47:06 JST 2011]
[am@aofd30 processor=x86_64]
[~/teaching/Lowpass/Lanczos]
$ make
if [ ! -d ./obj ]; then \
mkdir -p ./obj ; \
fi
gfortran -c -c -o obj/lanczos_cosine.o lanczos_cosine.f
gfortran -c -c -o obj/lanczos.o lanczos.f90
gfortran -c -c -o obj/print_header.o print_header.f90
gfortran -c -c -o obj/print_ioinfo.o print_ioinfo.f90
gfortran -o lanczos ./obj/lanczos_cosine.o ./obj/lanczos.o ./obj/print_header.o ./obj/print_ioinfo.o
[~/teaching/Lowpass/Lanczos]
$ cat para.txt
¶
infle="input/sf_eddy_data.txt"
ofle="output/lp_sf_eddy_data.txt"
dt=10.
prd1=120.
&end
[~/teaching/Lowpass/Lanczos]
$ chmod u+x *.sh
[~/teaching/Lowpass/Lanczos]
$ run.sh
Input: input/sf_eddy_data.txt
Output: output/lp_sf_eddy_data.txt
Number of the input data, N= 100
[~/teaching/Lowpass/Lanczos]
$ cd plot
[~/teaching/Lowpass/Lanczos/plot]
$ chmod u+x *.sh
[~/teaching/Lowpass/Lanczos/plot]
$ ls ../output/
lp_sf_eddy_data.txt
[~/teaching/Lowpass/Lanczos/plot]
$ lg.sh ../output/lp_sf_eddy_data.txt
Bash script ./lg.sh starts.
Input : ../output/lp_sf_eddy_data.txt
Output : ../ps/lp_sf_eddy_data.ps
pstext: Record 6 is incomplete (skipped)
Done ./lg.sh.sh
ローパスフィルタのサブルーチン(lancos.f)
subroutine lanczos_cosine
$(x,y,dt,MX,NDAT,PRD1,ntrim,monitor_on_off)
入力
x: 入力データ
dt: sampling interval
MX: 入力データの配列サイズ
NDAT: 入力データの個数
PRD1: cut-off period
monitor_on_off: monitor_on_offが1ならば、フィルタ係数(weighting.txt)と周波数特性(frf.txt)を出力する
出力
y: 出力データ
ntrim: 両端でカットされるデータの個数
C***********************************************************************
C .f
C
C Lanczos-Cosine filter
C
C NOTE
C Lanczos-Cosine filter (LCF) passes ~ 50% of the power of the signal
C whose period equals the cut-off period. See Fig. 5.10.12 on p. 537
C of Emery and Thomson (1998). If you want to remove the signal whose
C period is longer than T1, cut-off period must be 2*T1. Note that
C the longer the cut-off period is, the more data points are lost from
C each end of the time series.
C
C References
C Filter design
C Emery, W. J. and R. E. Thompson (1998): Data analysis methods in
C phyiscal oceanpgraphy, 2nd Ed., Elsevier, 638pp.
C
C Filter response factor
Cj 岩田彰・北村正・横田康成(1995):ディジタル信号処理. コロナ社, 163pp.
cj (Page 62-63)
C
C Brief introduction to the digital filter used in physical oceanogrphy
Cj 花輪公雄・三寺史夫(1985): 海洋資料における日平均値の作成について-日
Cj 平均潮位を扱う際の留意点-.沿岸海洋研究ノート, Vol. 23, No. 1, 79-
Cj 87.
C Hanawa and Mitsudera (Bulletin on Coastal Oceanography, 1985)
C***********************************************************************
subroutine lanczos_cosine
$(x,y,dt,MX,NDAT,PRD1,ntrim,monitor_on_off)
C------------------------- DECLARATION ---------------------------------
C implicit none
C implicit double precision(a-h,o-z)
C------------------------- INCLUDE FILE --------------------------------
C include 'comblk.h'
C-------------------------- PARAMETERS ---------------------------------
parameter(MW=10000)
PARAMETER(PI=3.14159265358979323846E0)
C------------------------- COMMON VARIABLES ----------------------------
C common //
C----------------------------ARGUMENTS ---------------------------------
C double precision
C Input
integer MX,NDAT, monitor_on_off
dimension x(MX)
real PRD1, dt
C PRD1: cut-off period
C dt: sampling interval
C Output
dimension y(MX)
integer ntrim
C character
C------------------------- LOCAL VARIABLES -----------------------------
C Weighting
dimension w(-MW:MW)
C Filter response factor
dimension R(MW*3),T(MW*3)
C Print date, hostname, cwd, and username
character dirinfo*144, userinfo*72, hostinfo*72
INTEGER date_time(8)
CHARACTER REAL_CLOCK(3)*12
C-----------------------------------------------------------------------
C---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
C Set number of the terms in weighting
ntrim=int(prd1/dt/2.*1.5) !number of distnct filter coefficients
if(ntrim.gt.MW)goto 8000 ! Error handling
if(ntrim*2.gt.NDAT)goto 8010 ! Error handling
rm=float(ntrim)
C Set weighting
* See Emery and Thomson (1998), pp.519-539.
omega_c=2.0*pi/(prd1)
omega_M=(rm-1.0)/rm
omega_N=pi/dt ! Nyquist frequency
w(0)=omega_c/omega_N
do k=1,ntrim
rk=float(k)
C sigma-factors to avoid overshoot ripples in filter response factor.
C See Fig. 5.10.12 on p. 537 of Emery and Thomson (1998).
sigma=sin(pi*rk/rm)/(pi*rk/rm)
CD
C sigma=1 !which means the cosine filter
CD
w(k)=(omega_c/omega_N)*sin(pi*rk*omega_c/omega_N)*sigma
$ /(pi*rk*omega_c/omega_N)
w(-k)=w(k)
end do !k
C Clear output variable
do i=1,NDAT
y(i)=0.
end do !i
C Apply filter
do i=ntrim+1,NDAT-ntrim
do k=-ntrim,ntrim
y(i)=y(i)+w(k)*x(i+k)
end do !j
end do !i
C
C monitor_on_off .eq.1 ==> Print filter specification
C
if(monitor_on_off.ne.1)then
return
end if
C Weighting
open(99,file='weighting.txt')
iunit=99
call print_header(99)
write(iunit,'(A,1x,i10)')'# ntrim = ',ntrim
write(iunit,'(A,1x,i10)')'# NDAT = ',NDAT
wsum=0.
do k=-ntrim,ntrim
wsum=wsum+w(k)
write(99,'(i10,1x,e12.4)')k,w(k)
end do
write(99,'(A,f12.4)')'# sum = ',wsum
close(99)
C Filter response factor (FRF)
nout=int(PRD1*5.0)/dt
istart=ntrim/4+1
do i=istart, nout
T(i)=dt*float(i)
R(i)=0.
sum=0.
do k=1,ntrim
rk=float(k)
omega=2.0*pi/T(i)
sum=sum+w(k)*cos(omega*rk*dt)
end do !k
R(i)=w(0)+2.0*sum
end do !i
open(99,file='frf.txt')
iunit=99
call print_header(99)
write(99,'(A,e14.6)')'# dt = ',dt
write(99,'(A)')'# T frf'
do i=istart,nout
write(99,'(f10.5, e14.6)')T(i), R(i)
end do !i
close(99)
return
C---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
8000 write(*,*)
write(*,'(A)')'Error: in subroutine lanczos_cosine'
write(*,'(A)')'ntrim violates dimension limit.'
write(*,'(A)')'ntrim must be smaller than or equal MW.'
write(*,*)
stop
8010 write(*,*)
write(*,'(A)')'Error: in subroutine lanczos_cosine'
write(*,'(A)')'ntrim is larger than number of the input data.'
write(*,'(A)')'ntrim must be much smaller than NDAT.'
write(*,'(A,1x,i10)')'ntrim = ',ntrim
write(*,'(A,1x,i10)')'NDAT = ',NDAT
write(*,*)
stop
end
lancos.f90 (上記のlancos.fをFortran 90に書き換えたもの)
!***********************************************************************
!
! Lanczos-Cosine filter
!
! NOTE
! Lanczos-Cosine filter (LCF) passes ~ 50% of the power of the signal
! whose period equals the cut-off period. See Fig. 5.10.12 on p. 537
! of Emery and Thomson (1998). If you want to remove the signal whose
! period is longer than T1, cut-off period must be 2*T1. Note that
! the longer the cut-off period is, the more data points are lost from
! each end of the time series.
!
! References
! See below.
!***********************************************************************
subroutine lancos(x,y,dt,MX,NDAT,PRD1,ntrim,monitor_on_off)
parameter(MW=10000)
PARAMETER(PI=3.14159265358979323846E0)
integer,intent(in):: MX,NDAT, monitor_on_off
real,intent(in):: x(MX)
real,intent(in)::PRD1, dt
! PRD1: cut-off period
! dt: sampling interval
real,intent(inout):: y(MX)
integer,intent(inout)::ntrim
dimension w(-MW:MW) ! Weighting
dimension R(MW*3),T(MW*3) ! Filter response factor
!-----------------------------------------------------------------------
! Set number of the terms in weighting
ntrim=int(prd1/dt/2.*1.5) !number of distnct filter coefficients
if(ntrim.gt.MW)goto 8000 ! Error handling
if(ntrim*2.gt.NDAT)goto 8010 ! Error handling
rm=float(ntrim)
! Set weighting
! See Emery and Thomson (1998), pp.519-539.
omega_c=2.0*pi/(prd1)
omega_M=(rm-1.0)/rm
omega_N=pi/dt ! Nyquist frequency
w(0)=omega_c/omega_N
do k=1,ntrim
rk=float(k)
! sigma-factors to avoid overshoot ripples in filter response factor.
! See Fig. 5.10.12 on p. 537 of Emery and Thomson (1998).
sigma=sin(pi*rk/rm)/(pi*rk/rm)
!D
! sigma=1 !which means the cosine filter
!D
w(k)=(omega_c/omega_N)*sin(pi*rk*omega_c/omega_N)*sigma &
& /(pi*rk*omega_c/omega_N)
w(-k)=w(k)
end do !k
! Clear output variable
do i=1,NDAT
y(i)=0.
end do !i
! Apply filter
do i=ntrim+1,NDAT-ntrim
do k=-ntrim,ntrim
y(i)=y(i)+w(k)*x(i+k)
end do !j
end do !i
!
! monitor_on_off .eq.1 ==> Print filter specification
!
if(monitor_on_off.ne.1)then
return
end if
! Weighting
open(99,file='weighting.txt')
iunit=99
write(iunit,'(A,1x,i10)')'# ntrim = ',ntrim
write(iunit,'(A,1x,i10)')'# NDAT = ',NDAT
wsum=0.
do k=-ntrim,ntrim
wsum=wsum+w(k)
write(99,'(i10,1x,e12.4)')k,w(k)
end do
write(99,'(A,f12.4)')'# sum = ',wsum
close(99)
! Filter response factor (FRF)
nout=int(PRD1*5.0)/dt
istart=ntrim/4+1
do i=istart, nout
T(i)=dt*float(i)
R(i)=0.
sum=0.
do k=1,ntrim
rk=float(k)
omega=2.0*pi/T(i)
sum=sum+w(k)*cos(omega*rk*dt)
end do !k
R(i)=w(0)+2.0*sum
end do !i
open(99,file='frf.txt')
write(99,'(A,e14.6)')'# dt = ',dt
write(99,'(A)')'# T frf'
do i=istart,nout
write(99,'(f10.5, e14.6)')T(i), R(i)
end do !i
close(99)
return
!---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
8000 write(*,*)
write(*,'(A)')'Error: in subroutine lanczos_cosine'
write(*,'(A)')'ntrim violates dimension limit.'
write(*,'(A)')'ntrim must be smaller than or equal MW.'
write(*,*)
stop
8010 write(*,*)
write(*,'(A)')'Error: in subroutine lanczos_cosine'
write(*,'(A)')'ntrim is larger than number of the input data.'
write(*,'(A)')'ntrim must be much smaller than NDAT.'
write(*,'(A,1x,i10)')'ntrim = ',ntrim
write(*,'(A,1x,i10)')'NDAT = ',NDAT
write(*,*)
stop
end subroutine lancos
! References
! Filter design
! Emery, W. J. and R. E. Thompson (1998): Data analysis methods in
! phyiscal oceanpgraphy, 2nd Ed., Elsevier, 638pp.
!
! Filter response factor
!j 岩田彰・北村正・横田康成(1995):ディジタル信号処理. コロナ社, 163pp.
!j (Page 62-63)
!
! Brief introduction to the digital filter used in physical oceanogrphy
!j 花輪公雄・三寺史夫(1985): 海洋資料における日平均値の作成について-日
!j 平均潮位を扱う際の留意点-.沿岸海洋研究ノート, Vol. 23, No. 1, 79-
!j 87.
! Hanawa and Mitsudera (Bulletin on Coastal Oceanography, 1985)