#############################################################################
# scalogram関数
# ・自動的に最高レベルを設定する
# ・レベル間分割数はユーザが設定する
# ・擬似周波数はlfの値
# ・ウェーブレット変換結果の絶対値の分布をスケーログラムという
#############################################################################

scalogram <-

function (x,		#入力信号
	   a,		#スケール

	  noctave = log2(nextn(length(x))),
			#最高レベルJ
			#nextn(n)は、正整数n以上の一番小さい２のべき乗を返す
			#length(x)は、配列xの長さ（要素の数）を返す
												 #log2 なので整数値を返す
		
############################################################################
#
	 nvoice = 24,	#レベル間分割数（ユーザが変更可能）
#
############################################################################
		  		  
	xlab = 'Time', ylab = 'Frequency [Hz]', 
	...)

{
    if (any(is.na(a)))			#scalogram(x,NA) とすると実行される
        a <- cwt(x, noctave, nvoice, plot = FALSE)	#上で設定した解像度でCWTを計算
 
    	T = deltat(x) * (length(x) - 1)				
		#（ サンプリング時間区間）＝（データ数−１）×（時系列データ間の時間間隔）

    	f = pretty(log2(nextn(length(x))) - 1 - (0:(nvoice * noctave))/nvoice)
		#[（最高レベル）−１−（レベル）] に対して数列を生成
		
	lt = pretty(time(x))
		#（信号xをサンプリングした時刻）に対して数列を生成
		
    	lf = (1/T) * 2^pretty( rev(log2(nextn(length(x))) - 1 - (0:(nvoice * noctave))/nvoice) )
		#rev関数で 数列[（最高レベル）−１−（レベル）]を大きいほうから並べなおす
		#２のべき乗を適当な間隔でつくる（縦軸の値をとるため）
		#（擬似サンプリング周波数）をつくる
		
    	f = pretty(rev(log2(nextn(length(x))) - 1 - (0:(nvoice * noctave))/nvoice))
		#擬似周波数をつくる。縦軸は下のほうが大きい値になる
				
    	bt = 1/diff(range(time(x)))
		#range関数で（信号xをサンプリングした時刻）の最小値と最大値を返す
		#diff関数で最大値と最小値との差分をとる
		#bt = サンプリング時間の逆数 = 実サンプリング周波数
				
    	at = -bt * min(time(x))
		#（−実サンプリング周波数）×（最初のサンプリング時刻）
		
    	bf = -1/diff(range(rev(log2(nextn(length(x))) - 1 - (0:(nvoice * noctave))/nvoice)))
		#[（最高レベル）−１−（レベル）]を大きいほうから並べた数列
		#最高レベルと最低レベルとの差をとって、逆数にする
		#マイナスがつくのは便宜上？
		
    	af = -bf * max(rev(log2(nextn(length(x))) - 1 - (0:(nvoice * noctave))/nvoice))
		#（レベル）を分割した数列

    filled.contour( Mod(a),  #複素数のModulus（絶対値）を返す
			xlab = xlab, ylab = ylab,	
			plot.axes = { axis(1, at + bt * lt, lab = lt)
					#横軸：実は引き算になっている
					axis(2, af + bf * f, lab = format(lf, dig = 1))
					#縦軸
    					}, ...)

    invisible(list(t = time(x), f = lf, cwt = a))
		#三次元配列（時刻、擬似周波数、展開係数）をコンソールに表示しないようにする
}







#############################################################################
# 入力信号
#   t : 時間
#		適当にseq関数で設定する
#
#	x : 入力信号
#		波形は数式ではなくても、測定データを読み込めばよい
#
############################################################################
library(Rwave)


#### 信号の生成ここから（ユーザが変更）#####################################

t = seq(0,		#信号の開始時刻
	1,		#信号の終了時刻
	len=512	#サンプリング周波数とあわせる
	)

x = 2 * sin(2*pi*16*t) * exp(-(t -.25)^2/.001)		#波形の生成
x = x + sin(2*pi*64*t) * exp(-(t -.75)^2/.001)

#### 信号の生成ここまで ####################################################


x = ts(x, 	#サンプリングデータ

############################################################################
#
	 deltat = 1/512	#サンプリング周波数。lenとあわせる
#				 ２のべき乗の数（2,4,8,16,...）でなければならない
#
############################################################################
	   
	   ) 
		
scalogram(x,NA)		#サンプリングデータに対してスケーログラムを描画