DCT結果の見方

「音の周波数測定（マイクから入力）」(このリンクは夜間は接続できません。) を作って見ました。これは、Silverlightアプリケーションにすることで、プログラムのダウンロードを意識せずに利用できるサンプルにしようと始めました。しかし、始めてみるとセキュリティやプライバシー保護のためらしく多くの制約があることを知りました。（制約はわかっても、何が（誰が）保護されるのかはわかりませんが）

この制約の中で、音声データがNellyMoser形式で取得されることについて考えて見ます。
一般的には、音声の圧縮には DCT（離散コサイン変換）が使われるようです。この結果からパワーの大きい周波数が知れればデコードする必要はありません。

下記に記載した内容を調べて、左図のような表示をする Silverlightアプリケーションを作って見ました。
（Cと刻んである音叉を打撃したところ）
「マイク入力のスペクトル」に置きましたが、このリンクは夜間は接続できません。

（古い話ですが、「 NSV Split(NSV動画ファイル分割/AVI変換/エラーFIX)」を作った時に、分割位置の確認用にVP3のキーフレームのデコードを書きました。ドット位置がパワー順に整列することを利用していたように記憶しています。これは、2次元のIDCTだったわけです。今回は、一次元のIDCTを書いたことになりますが、この入力から周波数分布がわかるのかと言う話です。）

１．FFTと比べてみる

とにかく DCT をして、FFT と比較してみます。「音について」「周波数の分解能は?」で使ったデータをDCTに掛けてみます。（そのスクリプトは、「R言語」「サンプル」にあります。）

100,300,800Hzのサイン波を1000Hzで1秒間サンプリングした状態をデータにしています。図は、上がFFT、下がDCTの算出値をプロットしたものです。上の図は、「音について」「周波数の分解能は?」にあるものと同じです。
データは1000個の実数です。
FFTでは、1000個の複素数が算出されます。
DCTでは、1000個の実数が算出されます。
FFTでは、サンプリングレートの半分までが得られ、算出値の前半だけが周波数（Hz)を表すのに使われます。
DCTの方は、算出値が何か、データの順序が何を示すのか良くわかりませんが、周波数成分を知ることはできそうです。
(R言語の dct() は、dttパッケージをインストールして使いました。)

２．なぜ DCT で周波数成分がわかるのか

「FFTで任意区間の周波数だけを計算する」で FFT の計算を追ってみましたが、DCT も同じ性質を利用しています。

通常サンプリングデータは、「サンプリングデートとサンプリング時間」、あるいは「サンプリングレートとサンプル数」の条件を持っていますが、ここではサンプリングレートは考えません。取り扱うのはサンプリングデータ数だけです。ここで求める周期、周波数はサンプリングデータ数（分の時間）に対して1周期、2周期と言う意味です。

知りたい周波数の周期を思った係数列をサンプルに掛けて和を求めます。その周期性がサンプリング・データあれば、符号が調整されることで和の絶対値が大きくなります。
FFTの場合はパワーが算出されますが、DCTの説明では単にDCT係数と言うようです。

左図は、8サンプルの場合の周期性を持った係数の例です。

(1,1,1,1,-1,-1,-1,-1)
(1,1,-1,-1,1,1,-1,-1)

の、2つの数列で、周期性があると積和が大きくなることを試してみます。

左図のように、8サンプルで2周期のデータを2つ用意しました。この二つは、初期位相が異なります。

まず、赤丸でプロットした8個のサンプリングデータに、それぞれ(1,1,-1,-1,1,1,-1,-1)を掛けて和を求めてみます。

(1,1,-1,-1,1,1,-1,-1)を掛けることで、符号がそろって、和が大きな値になります。

サンプル	0.000	0.975	-0.434	-0.782	0.782	0.434	-0.975	0.000	積計
基底	1	1	-1	-1	1	1	-1	-1
積	0.000	0.975	0.434	0.782	0.782	0.434	0.975	0.000	4.382

次は、(1,1,1,1,-1,-1,-1,-1)を掛けて和をとった場合です。

サンプル	0.000	0.975	-0.434	-0.782	0.782	0.434	-0.975	0.000	積計
基底	1	1	1	1	-1	-1	-1	-1
積	0.000	0.975	-0.434	-0.782	-0.782	-0.434	0.975	0.000	-0.482

加算によって値が小さいものになります。
このことから、サンプリングされたデータには、サンプリング数に対して2周期の成分が多く、1周期の成分が少ないと判断できるということだと思います。

もし、位相がずれていれば、符号が上手く変えられないようにも思いますが大丈夫なようです。
下記は、位相をずらしたデータでの計算です。絶対値では大きな差が付きます。

サンプル	0.866	-0.680	-0.563	0.931	0.149	-0.997	0.295	0.866	積計
基底	1	1	-1	-1	1	1	-1	-1
積	0.866	-0.680	0.563	-0.931	0.149	-0.997	-0.295	-0.866	-2.191

サンプル	0.866	-0.680	-0.563	0.931	0.149	-0.997	0.295	0.866	積計
基底	1	1	1	1	-1	-1	-1	-1
積	0.866	-0.680	-0.563	0.931	-0.149	0.997	-0.295	-0.866	0.241

実際の DCT の変換行列を見ます。

DCT の式を左のように表すとき、

変換行列の各行を図にすると、

この図と下記のR言語のスクリプトの実行結果からして、3行目 i=3 が1周期のパワーを表しているものと考えます。
i=5 が2周期です。FFTの結果は、半分で折り返しで共役数が入り、半分だけが使われますが、この場合は、奇数だけを使うので有効範囲は同じになります。
値の大きさは、FFTの半分程度のおおきさです。（スクリプトの最後の「FFTでDCT」から）
また、符号が付いていますが、パワーは絶対値で大小を判断します。

> #基底の計算
> N<-8
> f<-function(ii,jj){
+ r<-cos(pi*ii*(jj+1/2)/N)
+ }
> z<-outer(0:(N-1),0:(N-1),f)
> #基底のバー表示
> mat <- matrix(c(1,2,3,4,5,6,7,8), 2, 4, byrow = TRUE)
> layout(mat)
> for(i in 1:8)
+ {
+ par(mar=c(2.5,2.5,2,1))
+ barplot(z[i,1:8],ylim=c(-1,1))
+ title(main=i-1)
+ }
>
> #DCTの算出
> s2<-c(0.000 ,0.975, -0.434, -0.782, 0.782, 0.434, -0.975, 0.000)
> for(i in (1:8))
+ {
+ r[i]<-round(sum(s2*z[i,1:8]),3)
+ }
> r
[1] 0.000 0.834 0.000 1.340 0.000 -3.382 0.000 -0.271
>
> #dttパケージのdct()
> round(dct(s2),3)
[1] 0.000 0.834 0.000 1.340 0.000 -3.382 0.000 -0.271
>
> #wikipediaのDCT-IIの記載の確認（FFTでDCT）
> ss2<-matrix(0,1,32) #32個全部ゼロ
> for(i in 1:8) ss2[i*2]<-s2[i] #2,4,6,8,10,12,14,16
> for(i in 1:8) ss2[32-i*2+2]<-s2[i]#32,30,28,26,24,22,20,18
> ff<-abs(fft(ss2))
> round(ff[1:8]/2,3)
[1] 0.000 0.834 0.000 1.340 0.000 3.382 0.000 0.271
>

３．なぜ DCT は圧縮になるのか

３．１．振幅データの圧縮との比較で

時間に対する振幅データ（PCMデータ）をそのまま圧縮する場合は、音声データの特徴を利用できない。たとえば16ビットPCM音声サンプルはすべての値を取るし、明瞭な偏りもありません。
これに対して、周波数に対するパワーのデータは、偏りがあり、人の聴覚性能も利用できます。

３．２．なぜFFTでないのか

音の分析の話では主にFFTが、音声圧縮ではDCTが説明されます。この理由はわかりません。
私にもわかるのは、以下のことですがパソコンで扱う限りでは特に使い分ける理由になるとは思いません。

実数演算で結果も実数
cos だけで良い。（sinを使わない）

４．MDCT（修正離散コサイン変換）

実際の音声圧縮コーデックで使われるのは、MDCTのようです。
DCTはサンプルを固定サイズのブロック（一定時間分のサンプル）に分割して周波数軸に変換します。これを複合するときにもブロック単位に計算されるので、複合されたブロック間の時間に対する振幅の波形の連続性は保存されていないことになります。(左図)
連続するブロック間で、それぞれに同じ周波数の成分があっても、位相の連続性はないことになります。

これを解決するのか軽減するのかはわかりませんが、MDCTでは入力を出力より長い時間のデータで行っています。(左図）
入力は同じデータが重複して使われます。MDCTには、重複するサンプル数や使用する窓関数などデコードするためにはいろいろなことを考えなければなさそうです。

しかし、DCTではなくMDCTだと言うことは、デコードせずに周波数成分を知ろうとすることには影響なさそうです。いずれでも、出力は同じ意味を持っています。

５．デコーダのどこのデータを使うか

作成したデコーダの処理を眺めて、どこのデータを使うのが良いか考えます。

５．１．テストデータ

デコーダの検証も兼ねて、内容の分かっているテスト用のPCMファイルを用意して置きます。
これを ffmpeg でエンコードして、デコーダの入力にします。

５．１．１．データの作成

下記、R言語のスクリプトを使いました。
125Hz の純音を、サンプリングレート：8000Hzで1秒サンプリングした状態です。（サンプル数を256の倍数になるように8192にしてあるので1.024秒）

source("WavIO_R.R")
#125Hzの音を8000Hzで1秒間サンプリングしたデータを作る
t <- seq(0,by=(2*pi*125)/8000,length.out=8192)
s <- sin(t)*32767
#WAVファイルに出力
wav16write("test.wav",s,hz=8000)
#振幅データの表示
windows(5,5)
par(mar=c(4,4,3,1))
layout(matrix(c(1,2),2,1))
plot(1:length(s),s,type="l",xlab="sample",ylab="振幅",main="テストデータ（全体）")
plot(1:512,s[1:512],type="l",xlab="sample",ylab="振幅",main="テストデータ（先頭の512サンプル）")

このスクリプトで作成した、test.wav を、ffmpeg で、FLV にエンコードします。これをデコードして、PCMに戻して確認しておきます。

C:\>ffmpeg.exe -i test.wav -acodec nellymoser nm.flv
C:\>ffmpeg -i nm.flv -acodec pcm_s16le pcm.wav

エンコード、デコードは、256サンプルを単位に行われるのでこの限りではほぼ再現性があり途中結果の確認もできそうです。
ただし、ffmpeg の nellymoser のエンコードは上手くいっていません。FLVの段階で終端の音量が大きくなっています。振幅が小さくなっていることを考慮してみる必要があります。
ffmpeg のデコード結果の終端の512サンプルを見ると左図のようになっています。

最後の256サンプル分は本来の振幅に復号されているようです。

最後のDCT結果と、それ以外のDCT結果の差から、DCTでのパワーの表現のスケールがわかりそうです。

ここで作成したFLVファイル(nm.flv)上には、256サンプル単位に、32個のパケットが格納されています。このパケットは、64バイトで表されています。この例では、すべて、同じスペクトル分布ですから、このパケットは同じ値になるはずで、確認したところ、最後の1つ以外は完全に同じ値でした。

５．２．デコーダの入力

256サンプルを 64バイトに圧縮するものらしく、デコーダの入力は、この64バイトのブロックになります。 Flashのマイク入力のデフォルトのサンプリングレート 8000Hz で考えると、 31.25ms分です。

下図は、前述のテストデータをFLVにしたものの音声データ部分です。ヘッダ16バイトとデータ64バイトの80バイトの繰り返しになっており、データ部分は１から31ブロックまで同じ値になっています。

1 (192) ===============
0 0 0 b3 8 0 0 41 0 0 0 0 0 0 0 52
c0 ff 0 14 18 63 8c 12 c6 18 63 8c 31 c6 b8 75
c0 df 75 dd 5d e3 ff d9 58 9a a5 95 59 99 a9 f7
f7 ff e7 ff ff ff e7 6f 1d f0 77 5d 77 d7 f8 7f
36 96 66 69 65 56 66 ea fd fd ff f9 ff ff ff f9
32 (2672) ===============
0 0 0 4c 8 0 0 41 0 3 e0 0 0 0 0 52
fe 86 63 31 e7 de 6b d0 3d f7 1c 84 e 42 d8 75
5d b1 77 5d f7 de 83 87 87 78 db b6 6d 5b 96 66
69 65 56 ff ff ff ff 1b fb 9b e4 14 72 f3 22 b6
4f b4 f ea a0 f 7a 78 78 78 78 38 33 73 66 66

５．３．DCT結果の復元

まず、圧縮を解除します。
エンコード時には、DCT結果の偏りを利用して、量子化、符号化することで、64バイトにしたわけですから、これを復元すればスペクトルが得られるはずです。
　実際には、MDICTであることから、64バイトから128x4の結果が得られます。前述のテストデータの場合、最後を除くと、この4つの128個の組は、すべて同じ値でした。
　本来、サンプル数の半分（サンプリングレートの半分）の結果が得られるものなので、最初の128個の組を使うだけで良さそうです。
256サンプル分の時間で0,1,...,127周期までの、スペクトルが算出されているわけです。サンプリングレートが、8000Hzなので、0,31.25,...,3968.75Hzに相当します。

　前述のテストデータと、250Hz、500Hzのサンプルで試してみます。左図は、128個の組の先頭の24個をプロットしたものです。

Amazon.co.jp ウィジェット

mikeo_410