ランチョス・コサインフィルタ

プログラムとサンプルファイル

• 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

&para

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)