作図法、軸、単位パワースペクトルを書くのに、FFTの結果から計算した値を、X軸（FFTの結果のインデクス）から立ち上がる線で描いていました。一般的なパワースペクトルの図は、X軸についていないぎざぎざのグラフのようです。
また、音高を変える変換を考えているのですが、途中で値の範囲がそろわなくて行き詰まってしまいました。
そこで、一般的な作図法、軸、単位を少しずつ調べながら進もうと思います。 データの種類 | データの名称 | 図の名前 | X軸 | Y軸 | 備考 | 窓関数 | 窓関数の効果 | 窓関数の値 | 配列の順序 | 振幅(0から1) | | 窓関数の値をFFTしたもの |
周波数
| パワースペクトル(dB) | | パディングしたPCMデータ | ゼロ・パディングの効果 | | | | | PCMデータ | 振幅音声波形 | | 時間(秒) | 単位なし dB dBm | | 音量 | | | dB dBm | | スペクトル | パワースペクトル | ペリオドグラム | 時間(秒) | 単位なしマグニチュード(dB)パワースペクトル(dB) | detrend taper
FFT 補正 | スペクトログラム | 時間(秒) | 周波数(Hz) | 色、3Dでパワースペクトルを表す | | | | | | | | | | | | |
１．計算が合わない理由- FFTやその逆変換の計算とサンプル数で割るタイミング
- サンプリングデータのレンジ
16ビットサンプルなら整数値。sin()などで波形データを作ると、-1≦ ≦１。これを加算すると、．．． - FFTの計算結果の先頭とN/2+1をどう扱うか
２．ペリオドグラムの計算 
R言語のspectrum()を参考にして、計算方法を知ろうと思います。デフォルトで、１）デトレンド（トレンド除去）、２）テーパー、３）FFTの計算、４）補正の順に計算されています。 左図は、以下のR言語の実行結果です。 > s <- c(sin(seq(pi/8,by=2*pi/8 *1,length.out=8)),matrix(0,8))
> spectrum(s) 8サンプルで1周期のサイン波をπ/8から1周期と、8個の0で、16サンプルのデータを作って、spectrum()で描画しました。 その次の図は、
spectrum(s*32767)
として、描画しました。 さらに、次の図は、下記スクリプトで描きました。 
#データを作る。(サイン波の8サンプル、8個のゼロ)
s <- c(sin(seq(pi/8,by=2*pi/8 *1,length.out=8)),matrix(0,8))
y <- s # yがPCMデータ
N <- length(y) # Nがサンプル数
### データから線形トレンドを除去する ###
t <- 1:N - (N + 1)/2 # N個の数列。N=4なら、-1.5,-0.5,0.5,1.5
a <- (12 * sum(y * t)) / (N * (N^2 - 1)) # 回帰直線の傾き a
b <- mean(y) # 回帰直線の切片 b
y1 <- y - a*t - b # 回帰直線分を減ずる
### テーパー ###
taper <- 0.1
y2 <- Taper(y1, taper)
### 算出されるスペクトルの数は、N/2。Nが奇数なら切捨て ###
Nout <- floor(N/2)
### 算出されるスペクトルの数分の周波数の数列を作る。全体が1周期 ###
freq <- seq.int(1/N, by=1/N, length.out = Nout) # ゼロからではない
### FFTを計算 ###
f <- fft(y2)
pg <- f * Conj(f) / N # パワースペクトルを計算
### FFTの算出値の有効部分の取り出しと補正 ###
pg <- pg[2:(Nout+1)]
spec <- Re(pg[1:Nout]) / (1 - (5/8) * taper * 2)
### 作図 ###
windows()
plot(freq, spec, type="l",axes=T, log="y") # Y軸は対数 
■ X軸
デフォルトでは、0から0.5と表示されています。サンプリング周波数を与えていないので、１Hzとして計算しています。FFTの計算結果の1/2が表示されます。厳密には、FFTの計算結果を0からインデクスすると、１からN/2が使われています。先頭が使われないことは理解できますが、N/2は、良く分かりません。通常は、N/2は該当周波数のスペクトルとしては扱われないと思います。
■ Y軸
データ値に依存したスケールで、単位は無いと思います。対数軸で描画されています。FFTの計算は、-1≦ ≦1の値をサンプルに掛けてサンプル数合算するので、サンプル数で割れば、元の値より大きくなることはありません。この値の2乗なので、図のようなスケールになります。
２．１．トレンド除去 デフォルトが、detrend=Tで、以下のような計算を最初にしている。
定常過程であると見なすので、回帰直線が水平になるような変換を行う。これを「トレンドを除去する」と言うらしい。回帰線は重心を通るので、これをサンプリングデータから減ずる。
回帰直線、y = ax + b は、最小二乗法で求める。 
ここでは、X軸は時間なので、ｘを ｔと読み替える。また、サンプリング間隔は等間隔で増加すれば他に制約はないので、総和がゼロで、増分が1の数列を使うことにする。したがって、その平均値は0となる。 サンプル数を、N、サンプルを 
時間経過を表す数列を 
とすると、 
ｔの増分が1なので、 
プログラミングをする都合から、まず、以下を計算する。 
これらから、 
トレンドを除去した結果を、 
とすると 
となる。 ２．２．テーパー cosine-bell、taperと言うようです。
左図は、100個の1を入力して、spec.taper()と、R言語で書いたTaper()関数の計算結果を比べたものです。 Taper <- function(x,taper){
N <- length(x)
p <- floor(N * taper)
by1 <- pi / p
a <- (1-cos(seq(by1/2, by=by1, length.out=p))) / 2
r <- x
r[1:p] <- r[1:p] * a
r[N:(N-p+1)] <- r[N:(N-p+1)] * a
return(r)
} 気が付いたのは、先頭に掛けられる値がゼロではないと言うことです。演算範囲の減衰割合は1/2より少なくなると思います。また、補正時には先頭のデータを切り捨てているのでこれも考慮が要ります。 ２．３．FFTFFTの計算結果をからパワースペクトルを計算する段階でサンプル数で割っています。
パワースペクトルは絶対値を取って2乗するのではなく、共役数を乗じているので複素数のままになっています。 ２．４．補正 最終的にパワースペクトルの値として使われるのは、(1 - (5/8) * taper * 2)で除した値になります。5/8なる数字の意味は、私には分かりません。 補正が必要な理由は、パワースペクトルの総計はPCMデータのエネルギーからきており、計算で変わることはないからだと思います。補正の対象は、テーパーによって減衰させた分と、表示している値がFFT算出値の部分であることが考えられます。
テーパーによる減衰分は、０．４５から０．５と言った値です。これは、係数が０から１ではなく、その内側が使われているため若干少なくなっています。
０から数えて、１からN/2を採用している影響は良く分かりません。先頭はPCMデータの平均を表すので除外することは理解できます。N/2は、元のPCMデータに、1、-1を交互に乗じたものの和なので、周波数成分を表すものではないように思います。 |
インデクス
