mikeo_410

 窓を掛ける

エネルギーが減る

  確かなのは窓を掛ければエネルギーが減ります。

  1. > sum(rep(1,1000)*hanning(1000))/1000
  2. [1] 0.4995
  3. > sum(rep(1,1000)*hamming(1000))/1000
  4. [1] 0.53954
  5. > sum(rep(1,1000)*blackman(1000))/1000
  6. [1] 0.41958

長い周期は化ける

  FFTの計算区間で半周期に満たないような波形は、窓を掛けるともっと短い周期になってしまうことが推測できます。同様に、1周期以上でも長いものは変形されて影響を受けます。
  左図の赤は、4000Hzでサンプリングした 1H の正弦波の40サンプル分です。
  青は、この40サンプルにハン窓を乗じた結果です。

  こうした周波数のバケは一部は、以下のように重ね合わせで改善されると言うことだと思います。左下の黒は、3回の重ね合わせの合算です。

重ね合わせ

  4000Hzでサンプリングした 440H の正弦波の最初の40サンプル分(10ms分)を、窓掛け、重ね合わせで延長してみます。
  左図は、120サンプル分を描いたものです。
  440Hzは40サンプルで4.4Hzで、40サンプルでは切りの良い値でないことが要点です。 
  120サンプルでは、13.2周期分になります。
  左中の図のように、最初の40サンプル分を繰り返すと、その接続箇所は連続しません。

  下の図は、最初の40サンプル分にハン窓を乗じたものを、20サンプルずつずらしながら加算して作成した120サンプルです。
  連続に見えますが、12周期弱になっています。

窓関数

窓関数の計算

  1. # ハン窓
  2. hann <- function(L){
  3.   w <- 0.5 - 0.5 * cos( 2*pi * c( 0:(L-1) ) / (L-1) )
  4.   return(w)
  5. }
  6. # ハミング窓
  7. hamming <- function(L){
  8.   0.54 - 0.46 * cos( 2*pi * c( 0:(L-1) ) / (L-1) )
  9. }
  10. # ブラックマン窓
  11. blackman <- function(L){
  12.   x <- c(0:(L-1))/(L-1)
  13.   0.42 - 0.5 * cos(2*pi*x) + 0.08 * cos(4*pi*x)
  14. }
  1. # 窓関数のプロット
  2. L <- 100 # 窓の長さ
  3. par(mar=c(3,3,2,1),mgp=c(1.5,0.5,0))
  4. plot(hann(L),type="l",ylim=c(0,1),col=2, xlab="インデクス",ylab="振幅")
  5. par(new=T)
  6. plot(hamming(L),type="l",ylim=c(0,1),col=3,axes=F,xlab="",ylab="")
  7. par(new=T)
  8. plot(blackman(L),type="l",ylim=c(0,1),col=4,axes=F,xlab="",ylab="")
  9. abline(v=50,col="gray")
  10. title(main="窓関数3種",cex.main=1)
  11. legend(35,y=0.3,c("ハン窓","ハミング窓","ブラックマン窓"),c(2,3,4),cex=0.8,bty="n")

窓のＦＦＴ

  1. # パワー(dB)
  2. # x:PCMデータ
  3. # pad:0を補ってFFTを計算する場合の0の個数
  4. db_power <- function(x,pad=0){
  5.   if(pad>0) p <- abs(fft(c(x,matrix(0,pad))))
  6.   else      p <- abs(fft(c(x)))
  7.   p <- p/max(p)
  8.   p <- ifelse(p<(10^(-16)),10^(-16),p)
  9.   db <- 20*log10(p)
  10.   return(db)
  11. }
  1. # 窓関数を比較します。(窓関数をFFT)
  2. L <- 128          # 計算する窓の長さ
  3. # 比較のために計算するFFTの条件
  4. m <- 8            # L * m サンプルでFFTを計算
  5. pad <- L * (m-1)  # 0を補う個数
  6. N <- L * (m/2) # FFTの計算結果の周波数の必要部分 
  7. #-----------------------------------------------------------------------
  8. # パワーを計算
  9. p_rect <- db_power(matrix(1,L),pad)
  10. p_hann <- db_power(hann(L), pad)
  11. p_hamming <- db_power(hamming(L), pad)
  12. p_blackman <- db_power(blackman(L), pad)
  13. # 1からN周期部分をプロット
  14. ymin <- min(p_hann[2:N],p_hamming[2:N],p_blackman[2:N])
  15. freq <- seq(0, by=L/(L*m), length.out=N)
  16. par(mar=c(3,3,2,1),mgp=c(1.5,0.5,0))
  17. plot(freq[2:N],p_rect[2:N],type="l",ylim=c(ymin,0),col="gray",
  18.      xlab="周波数(ナイキスト周波数が64)",ylab="振幅(dB)")
  19. par(new=T)
  20. plot(freq[2:N],p_hann[2:N],type="l",ylim=c(ymin,0),col=2,
  21.      axes=F,xlab="",ylab="")
  22. par(new=T)
  23. plot(freq[2:N],p_hamming[2:N],type="l",ylim=c(ymin,0),col=3,
  24.      axes=F,xlab="",ylab="")
  25. par(new=T)
  26. plot(freq[2:N],p_blackman[2:N],type="l",ylim=c(ymin,0),col=4,
  27.      axes=F,xlab="",ylab="")
  28. title(main=sprintf("窓関数(%dサンプル、%dbyte/FFT)",L,L*m)，ｃｅｘ．ｍａｉｎ＝1)
  29. legend(35,y=10,c("矩形窓","ハン窓","ハミング窓","ブラックマン窓"),
  30.                 c("gray",2,3,4),cex=0.7,bty="n")
  31. title(main="窓関数のＦＦＴ",cex.main=1)

  窓関数の数列そのものをＦＦＴした結果です。
    下図は、m=1 として、同じスクリプトで描いたものです。
  上の図の縦じまは、ゼロ詰によって生じているものです。

  また、矩形窓は描画されなくなっています。これは、値が常に-320dBになっているためです。-320dBと言う値は、小さすぎる値をこの値にしているためで、本来、検出される周期性がありません。

  こうした図を見かけましたが、窓には周波数のフィルタ機能があるのでしょうか。
  全周波数の含まれるサンプルがあれば、これに窓を掛けて見ればわかります。

  1. # テストデータの作成
  2. fi <- rep(complex(re=0,im=0),100)
  3. fi[seq(2,51,by=2)] <- rep(complex(re=0,im=1),25)
  4. y <- Re(fft(fi,inverse=T))
  5. y <- y/max(abs(y))
  6. par(mar=c(3,3,2,1),mgp=c(1.5,0.5,0))
  7. mf <- par(mfrow=c(2,1))
  8. plot(0:99,y,type="l",ylab="振幅",xlab="サンプル")
  9. title("テストデータ",cex.main=1)
  10. z <- abs(fft(y)[1:51])
  11. plot(0:50,z,type="l",col="blue",xlab="周波数",
  12.      ylab="パワー",xlim=c(0,50),ylim=c(0,3.3))
  13. title("テストデータの周波数分布",cex.main=1)
  14. par(mf)

  FFTの逆変換で、振幅波形を作って見ました。
  全周期にパワーを設定してもOKそうに見えましたが、窓を乗じてFFTを計算するとパワーが1周期と50周期に集中してしまい目的を達しませんでした。（下図左）
  そこで、奇数周期にだけパワーを設定してテストデータの振幅波形を作りました。（下図右）

全周期にパワーを設定すると 奇数周期にパワーを設定すると

  100サンプルで、窓を乗じて、１から49周期の奇数周期のパワーを見ます。

  結果は、左図のようになりました。
  まず、窓関数を乗じることは、周波数応じたフィルタの効果はないようです。

  いずれの場合も、窓を乗じたことで周波数の検出能力が落ち、奇偶によるパワーの差が潰れています。hamming()の場合は、いくらかは保存されているようです。
 

mikeo_410@hotmail.com