フィルタ処理（移動平均，自己相関係数による周期性の観察）
簡単なフィルタ処理
オリジナルの南方振動指数には、短周期変動が混在している。
ここでは、簡単なフィルタ処理を行い、比較的長周期の変動を抽出してみる。
今回扱うフィルタ処理は、「移動平均」である。
処理は2段階に分けて行う。
（１）121か月（10年分）移動平均の計算
はじめに、121か月（つまり10年）の移動平均を行う。
時系列データの順番をNで表す場合、
N=1～121の平均、N=2～122の平均、N=3～123・・・
と121個分のウィンドウを移動させながら、平均を求めていく。
この処理を行うことで10年以内の周期を取り除く。
※現代においてはあまり用いない。フーリエ変換などをベースにフィルタ処理を行うのが一般的
（２）SOI時系列から121か月（10年分）移動平均を引く
次に、オリジナルのSOI時系列から、121か月移動平均の成分を除去する。
この処理により、SOI時系列科10年以上のゆったりとした長周期変動を取り除く。
この時点でSOI時系列データは、10年以内の変動成分しかもっていない、としたい。
※121か月移動平均のデータは前後60か月分、つまり120か月分のデータを失っていることに注意。
（３）さらに13か月の移動平均を行う
次に13か月移動平均を施すことで、1年以内の変動成分（例えば季節変動）を除去する。
サンプルコードはここ
GMTで描くSOI時系列（オリジナルと処理後）
gmt begin SOItimeseries
gmt set FONT_HEADING 24p
gmt basemap -R1951/2023/-4/4 -JX15c/7c -Bxaf5+l"Date" -Byaf0.5+l"Southern Oscillation Index"
gmt plot SOItimeseries.txt -R1951/2023/-4/4 -W0.2p -Bxag5 -Byag0.5 -BWSne+gazure1 -lSOI
gmt plot SOItimeseries.txt -Sc0.1c -Gwhite -W0.2p,black
gmt plot SOIsmoothed.txt -W0.8p,red -lSOI_smoothed
gmt legend -DjBL+o0.3c -F+gwhite+pthicker --FONT_ANNOT_PRIMARY=10p,Helvetica-Bold
gmt end show
自己相関係数
まず相関係数を理解しよう。ベクトルがわかっていれば問題なし。
相関係数についてはここを参照
標準偏差と共分散がわかっていればどうってことはない。
変数XとYの相関係数は、XとYの共分散をそれぞれの標準偏差で割ればよい。
自己相関を求めるコード（ラグ相関）
このままでは同時相関しかも求められません。
c============================================
C This code reads sea level purresure dataset
C from NOAA & cal Southern Oscillation Index!
C K.Mizobata
C 2024.03.28
c============================================
C234567--------------------------------------
character*120 header
real*8 date(732),SOI1(732)! smoothed data
real*8 SOI2(732)! smoothed data
real*8 COR(160) ! for correlation coeff.
!------ input "smothing" data ------!
open(10,file="SOIsmoothed.txt")
!do k=1,732!
do k=1,576 ! until Jun 2014
read(10,*) date(k),SOI1(k)
SOI2(k)=SOI1(k) ! copy!!
enddo
close(10) !<-- caution: should be closed
!=============================
! START CAL AutoCorrel
open(20,file="AutoCor.txt")
do i=1,160 !ラグを制御
!=============================
!------ cal mean -----!
icount=0
ave1=0; ave2=0;
!do k=1,732 ! all data
do k=1,576 ! until Jun 2014
ave1=ave1+SOI1(k) ! 合計の計算
ave2=ave2+SOI2(k)
icount=icount+1
enddo
ave1=ave1/real(icount) ! 平均の計算
ave2=ave2/real(icount)
!------ cal mean -----!
var1=0; var2=0; covar=0
std1=0; std2=0;
!do k=1,732
do k=1,576
var1=var1+(SOI1(k)-ave1)**2
var2=var2+(SOI2(k)-ave2)**2
covar=covar+(SOI1(k)-ave1)*(SOI2(k)-ave2)
enddo
var1=var1/real(icount) ! 分散
std1=sqrt(var1) ! 標準偏差
var2=var2/real(icount)
std2=sqrt(var2)
covar=covar/real(icount) ! 共分散
COR(i)=0 ! ゼロリセット
COR(i)=covar/std1/std2
!write(*,*) ave1,ave2,var1,var2,std1,std2,covar
write(*,*) "COR(",i,")=",COR(i)
write(20,*) i,COR(i)
!=============================
! FINISH CAL AutoCorrel
enddo
close(20)
!=============================
stop !<---- we need to add it ALWAYS!
end
GMTで自己相関係数
gmt begin AuoCor
gmt set FONT_HEADING 24p
gmt basemap -R1/160/-0.5/1 -JX15c/7c -Bxaf5+l"Lag(month)" -Byaf0.5+l"Correlation coeff."
gmt plot AutoCor.txt -R1/160/-0.5/1 -W0.2p -Bxag5 -Byag0.5 -BWSne -lCorrelation
gmt plot AutoCor.txt -Sc0.1c -Gwhite -W0.2p,black
gmt legend -DjTR+o0.5c -F+gwhite+pthicker --FONT_ANNOT_PRIMARY=10p,Helvetica-Bold
gmt end show
