Wavelet and short-time Fourier transform.

Continuous Wavelet transform with R:

折角なんだからMATLAB 使えよって気もするけど、入手しやすいRで。

R での Wavelet Transform に関する情報は Google 先生に聞いてもよくわからなかったので、ちょっと簡単な例をためしてみる。(legendにバグあり)

require('Rwave')

fps = 15
nframes = 512
t = (1:nframes)/fps
seq01 = (1:nframes)/nframes

f = 1.0*sin(2*pi*t) + 0.1*rnorm(nframes) + 0.8*sin((0.6+seq01*0.8)*2*pi*t)
ts_f = ts(f, deltat=1/fps)

plot(ts_f, type='l')
spectrum(ts_f)

#freqwave(ts_f, NA, 8, color.palette=heat.colors)

a = cwt(ts_f, noctave=8, nvoice=48, plot=FALSE)
filled.contour(Mod(a), color.palette=terrain.colors, nlovels=32)

a = cgt(ts_f, nvoice=256, scale=64, plot=FALSE)
filled.contour(Mod(a), color.palette=terrain.colors, nlevels=32)

matplot(cbind(Re(gabor(512, 256, frequency=2/15, 64)), 0.004*cos(2*pi*t)), type='l')
legend(0, 0.006, col=c(1, 2), lty=c(1, 2), legend=c("gabor", "cos 15Hz"))

まずは適当に信号(ts_f)を作成。
plot してみる。

spectrum してみる。これは多分、中身はfftだと思う。調べてない。

cwt してみる ... これは　'Rwave' ライブラリの関数だ。'wavelets' は dwt, 'Rwave' は cwt とかぽい。

Short-time Fourier transform with R:

時 間窓に Gauss 窓を使った Gabor 変換だ。ちょっと用語には自信なし。Wavelet 変換としての Gabor と Short-time Fourier transform としての Gabor って使い分けかたがあるのかね。まあどっちでも良いが、ここでは STFT だ。

cgt がそれ、のはず。なんか自信なくなってきた、、、。

周波数 1HZ (0.133)と、周波数が徐々にあがっていく信号が見える。えらいお気楽だな。
大学の授業では　DWT を C かなんかで書かされた記憶があるが。
ていうか、まず先に discrete じゃなくて continuous じゃなかろうか。聞いてなかっただけかな。

いまいち使い方に自信がない、ていうか、ドキュメントがなっていない。
ちなみに使われている Gabor は、こんな感じ(のはず)。@1Hz (0.133)

legend の 15Hz は間違えた。正しくは1Hzだ。1Hz, 15fps, 512 frame という設定のつもりだった。
今回は、一定周波数の信号、緩い周波数変化のある信号を見たかったので、時間窓をちょっと広めにとってみた。cwt より時間分解能が下がっているはずだが、まあ欲しい情報にあわせてできるだけ広くとるのが個人的には好みだ。

感想：
cwt は非常にお気楽、設定項目が少ない。Morlet が自動的に使われる。変更はできなっぽい。今回の例は簡単だからいいが、ノイズの多い信号の cwt を睨むのは結構大変そうだ。「うわーノイズだらけやな」とか感心する用途にはいいが、図から隠れた信号を見極めるのは職人芸ぽい。
一方 cgt はより図的には直感的。どういう信号があるはずで、どういうノイズが乗っているのかわかっていれば、自分で適当にいい感じの時間窓を選んであげれば、より直感的な図が得られると思った。

補足：
freqwave ... コメントアウトされてるのは google 先生に教えてもらった関数だから。ググるとちょっといい pdf が見れるかもしれない。
Emacs Speaks Statistics ... 良い感じです。
quartz (作図デバイス) ... 複数 png とかを出力しようとすると、絵が上書きされます。字とかAnti Alias とかきれいなのに。がっかり。ていうか、pdf が普通にGoogle Docs とかに貼れればいいんだけどな。