Rで経済時系列コミュのヒルベルト変換

桧木さんと電話で話していて表題の話になりました。

全然別の道楽で必要だったのですが、話しているうちに時系列予測に使えば
おもしろいんじゃないの？と言う話になり、Rで実装してみました。

とりあえず関数化したので貼っておきます。

中身を見て貰えばわかると思いますが、データを一旦fftし、半分から左をi倍し
半分から右を-i倍し逆fftで元に戻す変換です。

テレビなんかもデジタル化しましたが、あれのコア部分は、この変換をfirで行うもの。

たぶん相場に使っても、おなじような恩恵があると思います。

どういうもので、どう使うと面白いのかは、徐々にやって行こうと思います。

hirbertt <- function(data)
{
len <- length(data)
fdata <- fft(data)/len
fdata[1:(len/2)] <- 1i*fdata[1:(len/2)]
fdata[(len/2+1):len] <- -1i*fdata[(len/2+1):len]
result <- Re(fft(fdata,inverse=TRUE))
return(result)
}

plot(hirbertt(sin((1:628)/100)))

コメント(8)

最初
全て
最新の40件

[1] mixiユーザー 02月03日 11:08

一応私は工学系で言えば電気屋さん系の流れの一人です。

電気屋さんが他の工学と感覚的に違う所は、他の工学は基本的に静的な
解析を基本としています。

つまり動きが無い状態で、どうなるのか？と言う話。

一方で電気屋さんは、乾電池なら別ですが、家で電球ひとつ付けるにしても
交流と言う動きのある上で電球を付けるわけで、すぐに動きの上でどうなる
のかが解析の対象になる所が違うわけです。

で、まぁ動きがあるので、すぐに複素数で考える習慣が基本的にあるのですが
相場なんかは複素数で考えたりしないわけです。

ヒルベルト変換は、入ってきたデータの位相を９０度ずらす変換です。

単純に９０度ずらしたいだけなら、遅延させれば良いわけですが、どのくらい
遅延させれば良いかは周波数に依存するわけです。

相場のデータなんて、そんな決まった周波数があるわけでも無いですから
９０度ずらせと言われても、困るわけですが、ヒルベルト変換は、それが出来る
変換で、先ほどの例では、sinをヒルベルト変換するとcosが出てくる所を
計算させました。

色々なデータを入れてみて90度ずれるとどうなるのか適当に遊んでみて
下さい。

どうして90度ずれると、うれしいのか？と言うと複素数で表現する場合
たとえばインピーダンスであれば

Z = R + iX

とRとXが直交した話にして扱うわけですが、この直交を作るのに９０度
位相をずらすわけです。

これを相場で言えばRが変化してるわけですが、見えないだけでX分もある
と考えるとX分も変化しているわけです。

diff(Nile)とhirbertt(diff(nile))は、R成分とX成分を各々表しています。

こういう複素座標で見れば、結局は位相が動いているのか振幅が動いて
いるのかわかる為、予想がたてやすくなる効果や、今時のテレビ内部で
やっているような信号処理が可能になり、結構あれこれ遊べそうです。

まぁ、こんな話を昨日したばかりでした。

[2] mixiユーザー 02月03日 12:50

とは言っても、そんな簡単にイイコト出来ないんじゃないの？と思いますか？

例えば以下を実行してみましょう。

dn <- diff(Nile)
hdn <- hirbertt(diff(Nile))
plot(dn,hdn)

ナイルを和分したものを横軸、それをヒルベルト変換したものを縦軸に
プロットしただけのものです。

けど、今まで見たことの無かった特徴が見えていませんか？

[3] mixiユーザー 02月03日 13:06

そんなの、どれでも同じだろ？と思いますか？

data <- rnorm(100)
hdata<-hirbertt(data)
plot(data,hdata,type="line")

こうやってみるのとは、かなり違うと思いますけどね？

[4] mixiユーザー 02月03日 15:28

ここらで一旦放置します。

はーん。と思い当たる事がある方いますか？

[5] mixiユーザー 02月04日 15:31

日頃気にすることも無いと思うけど２点だけ。

１）直交してること

AとBが直交してるかどうか？は内積が0になるかどうかです。

A %*% B

を計算して0になれば直交。実際にやってみると結構残ります。
FFTって如何にも完璧そうですが、実態は結構ずさん。
まぁ元の大きさ|A|なんかよりは、全然小さいですけどね。

2)解析可能性

周波数軸で見た場合、A+iBの負の周波数成分がゼロになることです。

ヒルベルト変換は、Aに対して、こういったBを求める変換です。

[6] mixiユーザー 02月09日 23:50

少しだけヒントを書いておきますね。

どうせこいつ電気屋流れだから、座標の中心からの大きさと角度のベクトルを計算するんだろう？
とは思わないで下さいね。

どうして、こうするとダメかと言うと、こうすると中心からの距離はF分布になる為に、とても予測が
面倒になります。

では、どうするのか？と言えば、点から次の点へのベクトルと捉えてやれば、ざっくりと解決します。

で、必要があれば、この時系列を和分すれば良いでしょう。

つまり最初から和分しなくても、この時点で和分しても良いと思います。

この先は、このベクトル量を予測すれば良いわけで、やりかたは散々やって来たので割愛します。

以上で終わりとしますね。

[7] mixiユーザー 02月11日 22:15

以下は、オマケです。

前回のものは、与えるデータ数が偶数の場合のみしか正確に
計算できません。

どうしてかと言えば、奇数個の場合、どまんなかも右半分と
みなして、-j倍するからです。

本来、奇数個の場合は、どまんかはゼロなので、そうなるように
修正しました。

また、ちょっと真剣にヒルベルト変換をfirで実装するハメになりそうで
firで実装する場合のタップ係数のチューニングを始めています。

解析用に使うなら、このプログラムで充分だと思いますが、制御に
使う場合、まずはfirの次数は出来るだけ少なくしないと、遅れが
酷くなり使い物にならなくなる為です。

けど、まぁ次数を下げると当然精度も悪化するので、同じ次数で
如何に性能を出すかが問題になるわけです。

さらに相手がパワーの無い組み込みなので、FFTなんかを毎回
計算なんか、とても出来ない為です。

そちらのほうは、ブログのほうでやっています。

http://blogs.yahoo.co.jp/nozawa_onsen_potta

以下は修正版

#新型ヒルベルト変換
hirbertt <- function(data)
{
len <- length(data)
fdata <- fft(data)/len
fdata[1:(len/2)] <- 1i*fdata[1:(len/2)]
fdata[(len/2+1):len] <- -1i*fdata[(len/2+1):len]
if((len/2-round(len/2)) > 0){
fdata[(len/2+0.5)] <- 0+0i
}
result <- Re(fft(fdata,inverse=TRUE))
return(result)
}

[8] mixiユーザー 02月24日 14:49

おまけその２

前回からも結構改良しました。

fftの出力をDC成分（一番左）とナイキスト周波数（まんなかあたり）の部分も、きちんと計算させるようにしました。

データが奇数個（ナイキストなし）でも偶数個（ナイキストあり）でも処理出来るようにしてあります。

また戻り値ですが、単純にヒルベルト変換した値を返すものから、解析信号を返すように変更しました。

半分あたりを計算させるのにas.integerを使っていますが、Rってこうしないと全てベクトルなので
ループ変数が少数化しちゃうんです。

しばらく、なんだかわからない現象でしたw

hirbertt <- function(r)
{
len <- length(r)
lenby2 <- as.integer(len/2)
f <- fft(r)/len
f[1] <- 0
if(len %% 2 == 0){
f[lenby2+1] <- 0
}else{
f[lenby2+1] <- f[lenby2+1]*(1i)
}
for( l in 2:lenby2){
f[l] <- f[l]*1i
}
for( l in (lenby2+2):len){
f[l] <- f[l]*(-1i)
}
x <- Re(fft(f,inverse=TRUE))
x <- x*1i
h <- r+x
return(h)
}

検算します。

h<-hirbertt(1:4)
Re(h)%*%Im(h)
[,1]
[1,] 0

h<-hirbertt(1:5)
Re(h)%*%Im(h)
[,1]
[1,] -1.110223e-16

ばっちりですw

ログインすると、みんなのコメントがもっと見れるよ