前処理（窓関数、ゼロ・パディング）前処理の必要性は、まだ理解できていません。図を描いてみたので保存しておきます。 １．窓関数窓関数を掛けることは、情報を失うことだと思うので、必要性が分かるまでは置いておこうと思います。失う情報量は40%程度にもなるように思います。 窓は、サンプルの列から、FFTの計算区間を切り出すものとして説明されていますが、ここではフレーム化したデータに対して適用する、フレーム中のサンプル数と同じ個数の数列として扱っています。 ハン窓、ハミング窓、ブラックマン窓について作図して見ました。 窓を計算する関数： #---------------------------------------------------------------------
# ハン窓 # L:窓の長さ hann <- function(L){ w <- 0.5 - 0.5 * cos( 2*pi * c( 0:(L-1) ) / (L-1) ) return(w) } #--------------------------------------------------------------------- # ハミング窓 # L:窓の長さ hamming <- function(L){ w <- 0.54 - 0.46 * cos( 2*pi * c( 0:(L-1) ) / (L-1) ) return(w) } #--------------------------------------------------------------------- # ブラックマン窓 # L:窓の長さ blackman <- function(L){ x <- c(0:(L-1))/(L-1) w <- 0.42 - 0.5 * cos(2*pi*x) + 0.08 * cos(4*pi*x) return(w) } #--------------------------------------------------------------------- # パワー(dB) # x:PCMデータ # pad:0を補ってFFTを計算する場合の0の個数 db_power <- function(x,pad=0){ if(pad>0) p <- abs(fft(c(x,matrix(0,pad)))) else p <- abs(fft(c(x))) p <- p/max(p) p <- ifelse(p<(10^(-16)),10^(-16),p) db <- 20*log10(p) return(db) } #--------------------------------------------------------------------- # FFTの結果を中心がゼロになるように並べ替える。 # [1]を直流、後半を負の周波数と考える場合、左が負で、右に行くほど周波数 # が大きくなるように入れ替える。xは、長さが偶数。 fftswap <- function(x){ r <- x n <- length(x) n2 <- n / 2 r[(n2+1):n] <- x[1:n2] r[1:n2] <- x[(n2+1):n] return(r) }
100サンプルのデータに適用する場合を想定した、各窓関数を計算してみました。
#-----------------------------------------------------------------------
# 窓関数を比較します。 L <- 100 # 窓の長さ #----------------------------------------------------------------------- # 1から(L/2)Hz部分をプロット plot(hann(L),type="l",ylim=c(0,1),col=2, xlab="インデクス",ylab="振幅") par(new=T) plot(hamming(L),type="l",ylim=c(0,1),col=3, axes=F,xlab="",ylab="") par(new=T) plot(blackman(L),type="l",ylim=c(0,1),col=4, axes=F,xlab="",ylab="") abline(v=50,col="gray") title(main="窓関数(100サンプル用に計算)") legend(35,y=0.3,c("ハン窓","ハミング窓","ブラックマン窓"),c(2,3,4))
|
|
|
pcm <- wav16record() # デフォルトは11025Hzで2秒間の録音です。
wav16write("c.wav",pcm) # とりあえず、保存します。 windows(4,4) plot(pcm,type="l",xlab="時間",ylab="振幅") # こんな波形です。 # 定数を計算しておきます。------------------------------------- len <- 11025 * 1.48 # 1.48秒に相当するサンプル数。 len8 <- round(len / 8) # 比較用に1/8から算出します。 d <- 11025 /len # 周波数分解能(FFTの結果の間隔) paste("サンプル数:",len,"比較サンプル数：",len8,"分解能：",round(d,3)) #-------------------------------------------------------------- s1 <- pcm[4000:(4000+len-1)] # 1.4秒分のサンプルをs1に s2 <- s1[1:len8] # 比較用に1/8(175ms分)をs2に p1 <- db_power(s1) # パワーを計算 p2 <- db_power(s2, len-len8) # ゼロ詰めしたもののパワーを計算 # 作図(周波数) ------------------------------------------------- b <- round(460 / d) # 460Hzから描画 e <- round(600 / d) # 600Hzまで描画 freq <- seq( d * b, by=d, length.out=(e - b + 1) ) # 周波数の軸用 windows(4,4) plot(freq, p1[b:e], type="l", xlab="周波数(Hz)",ylab="パワー(dB)") par(new=T) plot(freq, p2[b:e], type="l",col=2, axes=F, xlab="", ylab="") title(main="ゼロ詰めの影響(1.4秒と175ms)") legend(535,y=-2,c("全サンプル","1/8サンプル"),c(1,2)) # 作図(詳細位置) ----------------------------------------------- max_index <- which.max(p1[2:(length(p1)/2)]) b <- max_index - 10 # 最大値の前後 e <- max_index + 10 # 10個を描画 freq <- seq( d * b, by=d, length.out=(e - b + 1) ) # 周波数の軸用 windows(4,4) plot(freq, p1[b:e], type="b", lwd=2, xlab="周波数(Hz)",ylab="パワー(dB)") par(new=T) plot(freq, p2[b:e], type="b",col=2, axes=F, xlab="", ylab="") title(main="ゼロ詰めの影響(1.4秒と175ms)詳細") legend(526,y=-1,c("全","1/8"),c(1,2)) # 分解能の比較 ----------------------------------------------- paste("ゼロ詰めなしで計算した周波数(Hz)：",max_index * d) # ゼロ詰めで8倍に p3 <- db_power(s1, length(s1)*7) max_index2 <- which.max(p3[2:(length(p3)/2)]) max_index2 d2 <- 11025 / length(p3) paste("ゼロ詰めして計算した周波数(Hz)：",max_index2 * d2) |
次に、ゼロ・パディングで分解能が上がるのかを見たいと思います。
前述のスクリプトの最後の部分を実行すると、以下の表示になります。
| > # 分解能の比較 ----------------------------------------------- > paste("ゼロ詰めなしで計算した周波数(Hz)：",max_index * d) [1] "ゼロ詰めなしで計算した周波数(Hz)： 522.972972972973" > # ゼロ詰めで8倍に > p3 <- db_power(s1, length(s1)*7) > max_index2 <- which.max(p3[2:(length(p3)/2)]) > max_index2 [1] 6191 > d2 <- 11025 / length(p3) > paste("ゼロ詰めして計算した周波数(Hz)：",max_index2 * d2) [1] "ゼロ詰めして計算した周波数(Hz)： 522.888513513514" > |
どちらが真実に近いのかは、結局分かりません。
そこで、純粋に計算の問題として、サンプルを作って、周波数の検出能力を比べて見ました。
| > # サンプリングレート1000Hzで、10秒間サンプリングする。サンプル数は10000。 > # 分解能は、0.1Hz > # サンプルは100, 100.1, 100.2, ... , 101Hz。 > s <- matrix(0,11,10000) > for(i in 1:11){ + # 1秒（1000サンプル）で、100Hzは100回転。0.1Hzは0.1回転。 + # それぞれ、2π/(1000/100)、2π/(1000/0.1)ずつ進む。 + s[i,] <- sin(seq(0, by=2*pi*(1/10 + (i-1)*1/10000), length.out=10000)) + } > # 全体から > d <- 1/10 # サンプリング時間が10秒なので分解能は0.1Hz > for(i in 1:11){ + p <- db_power(s[i,]) + k <- which.max(p[2:(length(p)/2)]) + print(paste("INDEX:",k, " /", k * d,"Hz")) + } [1] "INDEX: 1000 / 100 Hz" [1] "INDEX: 1001 / 100.1 Hz" [1] "INDEX: 1002 / 100.2 Hz" [1] "INDEX: 1003 / 100.3 Hz" [1] "INDEX: 1004 / 100.4 Hz" [1] "INDEX: 1005 / 100.5 Hz" [1] "INDEX: 1006 / 100.6 Hz" [1] "INDEX: 1007 / 100.7 Hz" [1] "INDEX: 1008 / 100.8 Hz" [1] "INDEX: 1009 / 100.9 Hz" [1] "INDEX: 1010 / 101 Hz" > # 2000サンプルから > d <- 1/2 # サンプリング時間が2秒なので分解能は0.5Hz > for(i in 1:11){ + p <- db_power(s[i,1:2000]) + k <- which.max(p[2:(length(p)/2)]) + print(paste("INDEX:",k, " /", k * d,"Hz")) + } [1] "INDEX: 200 / 100 Hz" [1] "INDEX: 200 / 100 Hz" [1] "INDEX: 200 / 100 Hz" [1] "INDEX: 201 / 100.5 Hz" [1] "INDEX: 201 / 100.5 Hz" [1] "INDEX: 201 / 100.5 Hz" [1] "INDEX: 201 / 100.5 Hz" [1] "INDEX: 201 / 100.5 Hz" [1] "INDEX: 202 / 101 Hz" [1] "INDEX: 202 / 101 Hz" [1] "INDEX: 202 / 101 Hz" > # 2000サンプルにゼロ詰めをして10,000サンプルに > d <- 1/10 # サンプリング時間が10秒なので分解能は0.1Hz > for(i in 1:11){ + p <- db_power(s[i,1:2000],8000) + k <- which.max(p[2:(length(p)/2)]) + print(paste("INDEX:",k, " /", k * d,"Hz")) + } [1] "INDEX: 1000 / 100 Hz" [1] "INDEX: 1001 / 100.1 Hz" [1] "INDEX: 1002 / 100.2 Hz" [1] "INDEX: 1003 / 100.3 Hz" [1] "INDEX: 1004 / 100.4 Hz" [1] "INDEX: 1005 / 100.5 Hz" [1] "INDEX: 1006 / 100.6 Hz" [1] "INDEX: 1007 / 100.7 Hz" [1] "INDEX: 1008 / 100.8 Hz" [1] "INDEX: 1009 / 100.9 Hz" [1] "INDEX: 1010 / 101 Hz" > |
サンプリングレートが1000Hz、サンプリング時間が10秒のデータを作ります。
サンプル数は、10,000個になります。
100Hzから101Hzまで、0.1Hz刻みのサイン波をサンプリングした状態を、s に作りました。
１．まず、全サンプルからパワーを計算した場合です。
サンプリング時間が10秒なので、分解能が0.1Hzとなり、正しく周波数を検出しています。
２．次は、2000サンプルの場合です。
サンプリング時間が2秒なので、分解能は、0.5Hzになり、3段階の値になりました。
３．同じ2000サンプルですが、ゼロを詰め10秒分のサンプル数と同じにします。
見かけ上の分解能は、0.1Hzです。
結果は、正しく周波数を検出しています。
| > # 800サンプルにゼロ詰めをして10,000サンプルに > d <- 1/10 # サンプリング時間が10秒なので分解能は0.1Hz > for(i in 1:11){ + p <- db_power(s[i,1:800],9200) + k <- which.max(p[2:(length(p)/2)]) + print(paste("INDEX:",k, " /", k * d,"Hz")) + }[1] "INDEX: 1000 / 100 Hz" [1] "INDEX: 1001 / 100.1 Hz" [1] "INDEX: 1002 / 100.2 Hz" [1] "INDEX: 1003 / 100.3 Hz" [1] "INDEX: 1004 / 100.4 Hz" [1] "INDEX: 1005 / 100.5 Hz" [1] "INDEX: 1006 / 100.6 Hz" [1] "INDEX: 1007 / 100.7 Hz" [1] "INDEX: 1008 / 100.8 Hz" [1] "INDEX: 1009 / 100.9 Hz" [1] "INDEX: 1010 / 101 Hz" > |
さらに欲を出して、800サンプル(800ms)分で計算して見ました。0.1Hzの精度があります。
これで、AudioFCプログラムを修正して見ます。
インデクス
