音声再生速度変換プログラム(2)

FFTの結果、元の振幅データを、サイン・コサインの和であらわすことができます。この式は、単なる時間の関数なので、合成に使われる正弦波は、終始一様（途中で振幅が変わったりしない）です。長いサンプリング時間を考えると不思議な気はしますが。
このことは、FFTの結果を操作する場合、合成に使われる正弦波の連続性（位相の連続性）を保つ必要があることを意味します。この位相の連続性は、FFTを行う単位で保つ必要があるのだと考えました。
したがって、周波数を2倍にした結果、時間が半分になったデータを、周波数が2倍のまま、元の時間に戻すには、合成しなおすと考えた方が良さそうです。

１．再度FFTのおさらい

FFTとその逆変換で、完全に相互に変換できることは、前に確かめました。しかし、合成を考えると、無音の区間はどう表現されているのか疑問に思いました。
そこで、16サンプルのサイン波の間に、ゼロを16個挟んで48サンプルを作って、考えます。

図を書くのに使ったR言語のスクリプトは、「R言語」「サンプル」に入れました。インデクスが1からなので見にくいですが。ゼロ始まりにする方法を知りません。

このサンプルを、FFTして、結果を逆変換して見ます。

 当然ながら復元します。

FFTの結果を調べます。

10,11,12,13
37,38,39,40
に、絶対値が大き目の値があります。この番号は、48サンプル中での周期になります。

まず、FFTの結果全部を使って合成して見ます。

i <- complex(im=1)
y <- matrix(0,48)
t <- 0
for(k in 1:48){
  for(j in 1:48){
    y[k] <- y[k] + Re(f0[j])*exp(i*(j-1)*t) + Im(f0[j])*exp(i*(j-1)*t)*i
  }
  t <- t + 2*pi/48
}

 FFTの結果で、絶対値が大きい8箇所で合成して見ました。左図の黒の線です。

10,11,12,13
37,38,39,40

次に、 9,41の9と41(9の折り返し位置)を加えて合成して見ました。左図の赤線です。

左図は、後半の合成を除いたものです、Y軸のスケールは半分になりますが、波形は同じです。

i <- complex(im=1)
y <- matrix(0,48)
t <- 0
for(k in 1:48){
  y[k] <-        Re(f0[10])*exp(i* 9*t) + Im(f0[10])*exp(i* 9*t)*i
  y[k] <- y[k] + Re(f0[11])*exp(i*10*t) + Im(f0[11])*exp(i*10*t)*i
  y[k] <- y[k] + Re(f0[12])*exp(i*11*t) + Im(f0[12])*exp(i*11*t)*i
  y[k] <- y[k] + Re(f0[13])*exp(i*12*t) + Im(f0[13])*exp(i*12*t)*i
  t <- t + 2*pi/48
}

左図は、三角関数に直して、三角関数で合成したものです。

FFTの結果の実部、虚部を、それぞれR,Iと書きます。インデクスはゼロから数えます。
f(t) = R[9]exp(I[9]t)+I[9]exp(I[9]t)
     + R[10]exp(i10t)+I[10]exp(I[10]t)i
     + R[11]exp(I[11]t)+I[11]exp(I[11]t)i
     + R[12]exp(I[12]t)+I[12]exp(I[12]t)i
     + R[36]exp(i36t)+I[36]exp(i36t)
     + R[37]exp(i37t)+I[37]exp(i37t)i
     + R[38]exp(i38t)+I[38]exp(i38t)i
     + R[39]exp(i39t)+I[39]exp(i39t)i

半分で折り返して、共役関係にある。
(半分で折り返しなので、exp(i36t)=exp(i(36-48)t)=exp(-I[12]t))

f(t) = R[9]exp(I[9]t)+I[9]exp(I[9]t)i
     + R[10]exp(i10t)+I[10]exp(I[10]t)i
     + R[11]exp(I[11]t)+I[11]exp(I[11]t)i
     + R[12]exp(I[12]t)+I[12]exp(I[12]t)i
     + R[12]exp(-I[12]t)-I[12]exp(-I[12]t)i
     + R[11]exp(-I[11]t)-I[11]exp(-I[11]t)i
     + R[10]exp(-I[10]t)-I[10]exp(-I[10]t)i
     + R[9]exp(-I[9]t)-I[9]exp(-I[9]t)i

exp(ix)+exp(-ix)=2cos(x)
exp(ix)-exp(-ix)=-2sin(x)i
を使って、三角関数に

f(t) = R[9][exp(I[9]t)+exp(-I[9]t)] + I[9][exp(I[9]t)-exp(-I[9]t)]i
     + R[10][exp(I[10]t)+exp(-I[10]t)] + I[10][exp(I[10]t)-exp(-I[10]t)]i
     + R[11][exp(I[11]t)+exp(-I[11]t)] + I[11][exp(I[11]t)-exp(-I[11]t)]i
     + R[12][exp(I[12]t)+exp(-I[12]t)] + I[12][exp(I[12]t)-exp(-I[12]t)]i
     = R[9]*2cos(9t) + I[9]*i*2sin(9t)*i
     + R[10]*2cos(10t) + I[10]*i*2sin(10t)*i
     + R[11]*2cos(11t) + I[11]*i*2sin(11t)*i
     + R[12]*2cos(12t) + I[12]*i*2sin(12t)*i
     = 2R[9]cos(9t) - 2I[9]sin(9t)
     + 2R[10]cos(10t) - 2I[10]sin(10t)
     + 2R[11]cos(11t) - 2I[11]sin(11t)
     + 2R[12]cos(12t) - 2I[12]sin(12t)

２．振幅の概観を知るには

振幅の概観として、振幅のグラフを上下から緩やかな曲線で挟み込んだものをイメージしていました。包絡、あるいは包絡線という言葉を思い出して調べて見ると、「与えられた曲線族と接線を共有する曲線」とありハズレてはいないようです。しかし、音声の話では、ケプストラムの話が書かれており、「LPCによるスペクトル包絡」、「LPCケプストラムによるスペクトル包絡」、「FFTケプストラムによるスペクトル包絡」と言ったものが使われるようです。必ずしも、「接する」ではないようですが、時間的に伸縮するための基準に使うのですから十分です。
ケプストラムの説明には、音源(声帯など)と、調音系(唇など)の特性を分けることを意味し、「フーリエ変換の対数をフーリエ変換」したものと書かれています。
とすると、ケプストラムからFFTを使って求まる包絡は、横が周波数で、PCMの概観ではないのかも知れませんが、試して見ます。

２．１．FFTケプストラム

前述の、48サンプルで、FFTケプストラムを書いて見ました。

黒：元の振幅波形
赤：FFTケプストラム

前のスクリプトで、f0にFFTの結果が入っているので

２．２．包絡の分離

 前のスクリプトのf2を、再度FFTします。f2も入力と同じ48個ですが、半分から折り返しなので、有効なのは24個です。このうち、前半の12個と後半の12個で、別々にFFTを計算します。
このうち、若い方のインデクスの12個から得られたのが、左図の赤線です。
調音系(唇など)の特性を表すのだそうです。
（人の声は、音源となる声帯などで周期の変化の少ない早い信号を作り、それを唇などでゆっくりと調音していると考えるものらしい。）

横軸は、周期なので、PCMの波形ではなく、FFTの結果のスペクトルに対応した包絡のようです。
左図の縦軸は、FFTの結果のスペクトルの絶対値です。

 後半の、12個からは、左図ができました。
音源(声帯など)の特性を示すのだそうです。

ただし、そもそものデータが声ではありませんが。

３．さらに、おさらい

３．１．スペクトル

これからは、FFTの結果（あるいは、その逆変換の結果）をスペクトルと言うことにしようと思います。
逆変換では、振幅が算出されるので、大元の振幅(PCM)データも、スペクトルと言うことにします。
実際には、周波数スペクトル、時間スペクトルと書くことにします。時間スペクトルをFFT（あるいは、その逆変換）をすると、周波数スペクトルが求まり、再度計算すると時間スペクトルが求まるものと考えています。
(Wikipediaで「周波数スペクトル」を引くと、周波数を横軸としたグラフを指すようです。ここで決めた「周波数スペクトル」は、「振幅スペクトル」というようです。「周波数スペクトル」は、まったく逆なのかもしれません。わたしのための混乱回避策で、横軸が、「周波数」か「時間」かで呼び分けたいと思います。)

スペクトルはエネルギーの多寡を示しているのだと思います。ただし、スペクトルは複素数です。
時間スペクトルは、スピーカーのコイルを流れる電流の向きと量を想像しています。電流はゼロでない限り、エネルギーがありますが、変化しないもの（直流）は音としては検出できません。時間スペクトルは、音量として認識されます。時間スペクトルは、実部の値を使います。単にFFTの逆変換を行えば、元の時間スペクトルに戻り、虚部はゼロになります。（そもそも、実部にPCMデータを、虚部にゼロを入れてデータを作ります。)

周波数スペクトルは、絶対値の2乗をパワースペクトルと言うので、強さを示すことになります。

３．２．時間

FFTの計算で扱っているパラメータは、振幅データの、１）振幅、2)データ個数です。サンプリングレートもサンプリング時間も使いません。したがって、FFTの扱う時間は実時間ではありません。
例として使っている「サンプリングデータのつもり」データから、三角関数の和を出しましたが、実際に数字を入れると以下のようになります。

f(t) = 0.46cos( 9t)-0.10sin( 9t)
      +0.32cos(10t)-0.06sin(10t)
      -0.30cos(11t)+0.06sin(11t)
      -0.36cos(12t)+0.08sin(12t)

ここで、tについて考えます。
tの前の数字は、FFTの計算結果の配列のインデクス(ゼロから数える)であり、周波数を示します。FFTの結果の[1]は、FFTの計算に使ったデータの範囲で、1周期の周波数スペクトルの値です。[2]は、2周期の周波数スペクトルの値です。
tも、「FFTの計算に使ったデータの範囲で」と言うことになります。もし、「FFTの計算に使ったデータの収集時間を1とすると」、と考え c とすれば、
t=2π・c
です。

３．３．FFTの計算結果の対称性

48サンプルでは、絵に書くのが大変なので、8サンプルの例で考えます。サンプル数の半分-1 の周期まで検出できるので、2周期、3周期のサイン波でデータを作ります。3周期のサイン波は、sin()の値を半分にして、振幅を1にします。最後に直流成分のつもりで、0.2を加算しておきます。

前述の、振幅を表す8サンプルをFFTします。
スクリプトのインデクスは1からなのですが、説明では0からとして書きます。
[2]、[6]に、絶対値0.5があり、周波数2が強さ1である。
[3]、[5]に、絶対値0.25があり、周波数3が強さ0.5である。
と、読めます。
[0]は、直流分として加算した0.2が反映されています。
「FFTで任意区間の周波数だけを計算する」に書いた通り、[0]は、全サンプルに、
exp(complex(re=0,im=0))
を、乗じて合計したものです。複素数のexp(0)は、1.0+0・iなので、単に合計を取った物になります。これをサンプル数8で割って、0.2になります。
[4]は、exp(complex(re=0,im=0))とexp(complex(im=-2*pi/8*4))を、それぞれ偶奇に掛けて合算します。exp(complex(im=-2*pi/8*4))=-1+0・iです。インデクスが奇数のデータの符号を反転して合計を取った物になります。

どのように利用できるか分かりませんが、[0]は、サンプリングデータ（PCMデータ）の合計です。件数で割れば、平均値を表しています。
インデクスがサンプル数/2の位置には、exp(0)とexp(-π・i)を交互に掛けて合算した値が入ります。0から数えて、偶数番のサンプルには１を、奇数番のサンプルには-1を掛けて合計します。
したがって、[0]+[サンプル数/2]は、偶数番号のサンプルの合計を2倍した値、[0]-[サンプル数/2]は、奇数番号のサンプルの合計を2倍した値になります。
インデクスがサンプル数/2の位置の値は、符号の変化の程度を表し、ゼロ交差の程度を表すことになります。偶数のグループと奇数のグループのサンプルのが、それぞれ、すべて同符号で、奇偶で符号が反対なら、絶対値が最大になり、[0]と同じ絶対値になります。それ以外は、正負の加算による相殺で、[0]より小さい絶対値になります。

サンプル数の1/2（8/2=4）のインデクスを挟んで、[1][7]、[2][6]、[3][5]とペアで周期1,2,3の強さを示し、ペアの値の関係は、共役(虚部の符号が逆)の関係になっています。

これを、FFTの逆変換すると、実部に、最初の振幅信号と同じ値が求まります。

 私の疑問は、逆変換ではなく、FFTの結果の共役数を求めてからFFTしても同じ結果になることです。FFTなのに、前述のようなサンプル数の1/2で折り返すようにはなりません。差は、虚部の数字にありそうです。
0から7の連番をデータとしてFFTを計算してみます。このとき、虚部には、同じ値(0,1,2,3)を詰めてみます。虚部が同じ値なら、「中央で折り返して共役」の関係が生じるようです。

サンプリングデータは、実数なので、虚部0の複素数と見なされますが、複素平面で考えると、実軸上の点となります。（スクリプトは、R言語の「サンプル」に入れました。）
それぞれの点に、FFTの計算と同様に指数関数を掛けて見ました。赤が、1周期目に相当する[1]の計算結果で、緑が、そのペアの[7]の計算結果です。（重なるポイントは、スクリプトの記述順で、最後に書いた色になっています。部分的にスクリプトを実行すれば確認できます。）
連番なので、絶対値が大きくなりながら、同じ角度で回っています。赤と緑の同番号は、共役関係にあることが分かります。ペアが共役になるのは、指数関数の適用の順序によっているようです。

これが、逆変換で戻るのは、再度、指数関数を適用することによるようです。FFTの逆変換は、FFTの結果に対して、共役数を求め、これをFFTします。
そこで、赤のデータに、[7]に適用した指数関数を、緑に[1]に適用した指数関数を、それぞれ掛けて見ます。
下図のように、元に戻りました。
本当は、FFTの計算では、赤の値が保存されるわけではなく、和が[1]の値となるわけなので、正しくありませんが。

 そこで、もっと単純な4サンプルの例を見て見ます。指数関数としては以下のものが使われます。指数関数については、常に、先頭は、exp(0)で1で、サンプル数/2の位置は、-1だと言うことです。

e0 = 1
e1 = exp(-2π/4・i)
e2 = exp(-2π/4 ・ 2・i) = -1
e3 = exp(-2π/4 ・ 3・i)

a0,a1,a2,a3のサンプリングデータ(実数)があるとします。

FFTの結果、F0,F1,F2,F3が得られます。FFTの計算は、次のように行われます。
F0 = ((a0・e0) + (a1・e0) + (a2・e0) + (a3・e0)) / 4
F1 = ((a0・e0) + (a1・e1) + (a2・e2) + (a3・e3)) / 4
F2 = ((a0・e0) + (a1・e2) + (a2・e0) + (a3・e2)) / 4
F3 = ((a0・e0) + (a1・e3) +(a2・e2) + (a3・e1)) / 4

F0,F1,F2,F3の共役数を、それぞれf0,f1,f2.f3とすると
FFTの逆変換の結果は、以下のI0,I1,I2,I3が得られます。
I0 = ((f0・e0) + (f1・e0) + (f2・e0) + (f3・e0))
I1 = ((f0・e0) + (f1・e1) + (f2・e2) + (f3・e3))
I2 = ((f0・e0) + (f1・e2) + (f2・e0) + (f3・e2))
I3 = ((f0・e0) + (f1・e3) + (f2・e2) + (f3・e1))

ここで、I1の計算に注目します。FFTの計算で中心から折り返しで共役関係になっていることが分かっているので置き換えます。(先頭とサンプル数/2の位置はそのまま)
I1 = ((F0・e0) + (F3・e1) + (F2・e2) + (F1・e3))
    = ((F0・e0) + (F1・e3) + (F2・e2) + (F3・e1))

指数関数の列は、e0,e3,e2,e1となっています。これは、F3のときに使われた指数関数の順番と同じです。

FFTの計算では、[0]と[サンプル数/2]は特別な意味がありますが、周波数のスペクトルを表す[1]から[サンプル数-1]は、[サンプル数/2+1]から[サンプル数-1]と共役の関係があるので、いずれか一方の値があれば十分です。値は複素数です。その逆変換は、元の振幅信号（PCMデータ）に戻るので、実数で、サンプル数全体が振幅を表します。
ただし、これは、FFTと、その逆変換が持っている性質ではなく、処理するデータが持っている性質のようです。虚部が全て0・iのデータをFFTすれば、共役関係を持った複素数データになり、この複素数データをFFTすれば、 虚部が全て0・iのデータに戻ると言うことのようです。

３．４．FFTの計算手順

FFTの計算では、時間経過 t は角度で表されます。最初のサンプルがt=0で、最後がt=(2π-2π/(サンプル数))です。サンプル数 N が偶数だけを考えますが、0から数えて、N/2 のサンプルは、t=πです。
回転に使用される指数関数exp(-t・i)はN個で、exp[0]=exp(-2π/N・0・i)、exp[1]=exp(-2π/N・1・i)、...、exp[N-1]=exp(-2π/N・(N-1)・i)です。
FFTの計算の結果は、N個の複素数の配列で、[0]は、全サンプルに、exp(-0・i)を乗じて合計したものになります。exp(-0・i)=1なので、単に合計したものになります。（R言語のFFTは合計ですが、サンプル数で割っているFFTプログラムもあります。）
[1]は、全サンプルに、順次、exp[0],exp[1]、...、exp[N-1]を掛けたものの和になります。
[2]は、全サンプルに、順次、exp[0],exp[2]、...、exp[N-2],...を掛けたものの和になります。
[N-1]は、全サンプルに、順次、exp[0],exp[N-1]、...を掛けたものの和になります。
FFTの計算の結果のインデクスを k とすると、exp[]のインデクス j は、j=0,j=j+k,...と使われます。ｊ≧N なったらj=j-Nとして巡回します。

左図は、0,1,2,3,4,5,6,7の8個の連番を、FFTと逆変換したものを表しています。
黒い数字が、連番データです。実数なので、虚部が0・iとして、実軸上の点となります。
緑がFFTの結果で、abc順に、FFTの結果の配列の順番になっています。b-h、c-g,d-fがペアで、サンプル中で、1,2,3周期のスペクトルを表します。ペアは、実部、虚部のがそれぞれの絶対値が、それぞれ同じ値で、虚部の符号が逆の共役の関係になっています。

この緑の値を逆変換すると赤の値になり、逆変換の結果の配列には、0,1,2,3,4,5,6,7が入ります。

 FFTの結果の配列のうち、先頭の[0]は、データの合計になります。左図では、サンプル数で割ってあるので、平均値になります。(左図の3.5+0・iの緑のa）

中央の[4]（サンプル数の半分）は、データのインデクスが偶数には、1を奇数には-1を乗じて合算したものになります。データが実数なら、虚部は0になります。左図の緑のeが該当しますが、これもサンプル数で割ってあります。
サンプル数の半分の位置の値は、振幅を示しているデータのゼロ交差数の度合いを表すことになると思います。全て同じ符号なら相殺して小さい値になります。

 ３．５．振幅

Windowsパソコンで普通に使われている、符号付16ビット/サンプルのデータを扱います。値は、16ビットの整数と同じで、-32768から32767範囲です。振幅の中心はゼロで、スピーカーを考えると、ゼロを中心に電流の向きが反転し、数値の絶対値に応じた電流が流れると言うことだと思います。ただ、マイクを考えると（カーボン、コンデンサ、...）、ゼロが中心になる仕組みが心配になりますが。

FFTの計算は、浮動小数点数(double)で行っています。最終的には、元の-32768から32767範囲の16ビット整数にしますが、このとき、入力のレンジと同じレンジになるように考えます。レンジが変わると、音量が変わります。
気になるのは、-32768で、どう解釈するかには流儀がある可能性があります。正負同じ範囲だと、-32767までで、-32768はハンパです。
FFTの計算は、指数関数(-1≦ ≦1)を、個々のデータ値に乗じて、個数分加算したものなので、FFTの計算結果の個々の値は、データのレンジの内側にあります。これは、入力データか、計算結果がデータの個数で割ってあることが前提です。（R言語のFFTは、割っていないので、fft1<-fft(samples)/length(samples)のように使います。）また、複素数なので、レンジの話は、絶対値が、と言うことになります。

考えやすくするために、32768割って、-1から1の範囲に正規化することが考えられます。（前述の理由で、-1≦ ＜1）
ここで例に使っている「サンプリングデータのつもり」も、正規化された状態です。
しかし、FFTの計算自体は、正規化不要です。FFT（データ個数で割る）したものを、逆変換すれば、入力データのレンジにかかわらず、完全に同じ値に戻ります。（これは、doubleの仮数部が6.5バイトあり、16ビット整数ならロス無く変換可能と言うことです。(C#でプログラミングします。)）
この件で気になるのは、log()です。このlog()は、複素数の絶対値に対して使用しているので、実数のlog()です。

正規化の有無で、符号が変わります。
しかし、試して見ると、算出されるケプストラムには、差がないようです。
ここでも、正規化は不要のようです。

ただし、ここにあげたケプストラムの図は、[2:48]として、先頭の値を捨てて描画してあります。
なぜ、先頭の[1]だけが、大きく違うのでしょうか。
その前に、そもそもなぜ[1]以外は、同じ値になるのかを考えて見ます。

符号付き16ビットPCM信号の最大値は32767なので、１から32767の対数を対数化すると、0から10.39718の範囲になります。
正規化した場合、絶対値の最小値は、1/32768です。-10.39721から0が対数化した場合の範囲になります。
ゼロですが、log(0)は、-Infと表示され負の無限大です。 FFTの結果に絶対値0が無い理由考えつかないので、ゼロは特別に処理する必要がありそうです。
正規化していない場合は、1未満は1に、正規化してある場合は1/32768未満は、1/32768に、置き換えれば良いように思います。特に、正規化した場合は負の対数は、完全なゼロではなくても、0＜ ＜1/32768以下の数字は、負の無限大に近い値になってしまいます。

正規化していないものを対数化した結果は、正規化したものを対数化した結果にlog(1/32768)を加算したものになることが分かりました。

サンプル中に2周期の計算を見て見ます。

 符号の混じったものの加算によって、レンジだけが影響すると言うことのようですが良く分かりません。

 [0]の差異が大きく異なる理由は、[0]の計算では、乗ずる値が全て1（exp(0)）なので、単にFFTの結果の値だけで決まります。この値が、対数化によって、正規化した場合は、負のみに、正規化しない場合は正のみになり、その和はまったく異なった値になります。

４．もうすこし

４．１．ケプストラムは実数か

FFTの逆変換の結果の値も、複素数になります。単にケプストラムと言うと、この実部の値を指すようです。
FFTの結果を、元の時間と振幅のデータに戻して見たときには、実数で疑問に感じませんでしたが。

４２．絶対値の効果 

前述48サンプルの例のデータから計算したFFTの結果のf0[10]を見てみます。
FFTの対象となる時間（サンプリング時間）で、9周期の周波数に対応する値です。
FFTの結果は、中央で折り返していて、f0[40]も、9周期の周波数を示します。

左図の赤がf0[10]です。緑は、f0｢40]で、f0[10]と共役の位置にあります。
f0[10]の数値を見てみます。

FFTの結果は複素数であり、大きさだけでなく位相の情報を含んでいます。
絶対値は、実数です。最大振幅を示す、実数に変換することになります。

 この後、対数をとって、実部に絶対値、虚部にゼロを入れた複素数が作られ、FFTの逆変換（ケプストラムの算出）が行われます。

４．３．対数をとる効果

FFTの計算の結果の各値は、個々のデータに－１から１の指数関数を乗じたものを、データ個数分加算したものなので、データ個数で割れば、元の振幅データのレンジの中の値になります。ただし、FFTの結果は複素数なので、直ちに比較できません。
前の方で、FFTの結果から、三角関数の和に変換しましたが、FFTの結果は、データ個数の半分で折り返しで共役関係を持つペアからできています。前述の例では、f0の10番目と40番目がペアでした。ここからは、インデクスをゼロから数えることにして10番目（スクリプトの表現のf0[10]）を[9]と書きます。
[39]を[9]の実部、虚部の値で表して、
f(t) = R[9][exp(I[9]t)+exp(-I[9]t)] + I[9][exp(I[9]t)-exp(-I[9]t)]i
       ：
     = 2R[9]cos(9t) + 2I[9]sin(9t)
       ：
と、なっています。この2R[9],2I[9]が、元の振幅データのレンジを超えないと言うことです。FFTの各値の実部、虚部のレンジは、半分になっています。もし、単一の周波数からなり、-1から1に正規化されている振幅データなら、FFTの結果は、その周波数の位置のペアに、それぞれ0.5が入ることになります。絶対値が0.5で、位相によって実部虚部に案分されます。

今回の振幅データは -1から1 の範囲なので、データ数で割ったFFTの結果の絶対値は0から1の範囲になります。したがって、この値の対数は、
の範囲になります。（もし、FFTの結果の絶対値がゼロだと無限大）

１．全て、負の数になる。
２．レンジが２から　と広くなる
３．FFTの計算結果の絶対値の増減は、同じ方向

４．４．FFTの逆変換の入力データ

左図の、下方の赤い折れ線が対数化した、逆変換の入力データ。
情報の黒い折れ線が、FFTの結果の絶対値を取ったものです。
ピンク、緑のマルは、FFTの結果の実部、虚部です。

左の縦の軸は、最終的な逆変換の入力データのものです。
右の縦の軸は、FFTの結果と絶対値のものです。目盛りに注意してください。2倍に強調した形になっています。

４．５．FFTの逆変換の効果

■周波数とスペクトルの関係のFFTの結果に対して、絶対値をとり、その対数をとって強調した結果を、逆変換します。もちろん、無操作で逆変換すれば、元の振幅信号に戻ります。
逆変換は、周波数とスペクトルの関係から、時間と振幅の大きさへの変換です。

ケプストラムの計算から想像できるのは、強調によって、高い周波数と低い周波数の振幅データは多く出力され、中間的な周波数成分に応じた周波数部分は少なく出力されると言うことです。

 左図の、黒はケプストラムです。縦軸の目盛りは、ケプストラムのものです。他の値は、演算方法が異なり、それぞれのスケールになっているので、形状だけ見てください。
逆変換(緑)は、何の演算もしないで、単に逆変換したものです。入力の振幅データが復元されています。
shift(赤)は、ケプストラムの結果に似るように、FFTの結果を、シフトとスケーリングで移動したものを逆変換したものです。（このデータは、周波数による強調の操作がされていません。）

■ケプストラムを再度FFTして、周波数のスペクトルを求めてみます。

左図は、絶対値の2乗によって、パワースペクトルとし、さらに0-1になるように、スケーリングしてあります。

 緑は、最初の振幅データを、FFT、FFTの逆変換したものなので、最初のFFTの結果と同じで、パワースペクトルを表しています。

 ケプストラムの演算を経た結果、黒のようにスペクトルが変化しています。
赤は、ケプストラムの演算と同じ領域になるように、データをシフトとスケーリングしたものの逆変換し、再度FFTしたものです。
図「log(abs(x))の効果」のような変換の影響が分かると思います。indexの24が中央なので、12のあたりが検出可能な周期の中心になります。このあたりを、境に長い周期性と早い周期性が分離されると言うことだともいます。
この図は、FFTの計算結果なので、半分で折り返しになっているので、24までの前半を見ます。

 ■ケプストラムを再度FFTしますが、今度は半分をゼロにしてから行います。時間-振幅の関係から、周波数-スペクトルの変換になルはずですが、対数振幅スペクトルに対応するケプストラム係数と言うようです。

左図の、緑は元になったPCMデータのFFTの計算結果の前半部で、絶対値です。
黒が、ケプストラムから、後半に0を詰めてFFTした結果で、絶対値です。
赤は、シフト・スケーリングしたものです。 

用語の理解

FFTの入力

非定常信号
アナログ信号からデジタル信号へ。標本化、量子化、符号化
A/D変換
PCM

FFTの出力

短時間スペクトル、スペクトログラム、フーリエスペクトル

短時間位相スペクトル

複素数なので、パワースペクトル、短時間パワースペクトル（短時間スペクトルの絶対値の2乗）、振幅スペクトル（パワースペクトルの平方根、短時間スペクトルの絶対値）

 ケプストラム
通常、実数のケプストラム。フーリエスペクトルの絶対値の対数を逆フーリエ変換したもの。

パワーケプストラムは、実数のケプストラム。単にケプストラム。パワーケプストラムから計算されるが、対数をとったら、2で割って逆フーリエ変換するので、同じもの。

複素ケプストラム

複素数 z = r・exp(i・t) の対数は、
log(z)  = log(r) + i・t

これを、逆フーリエ変換したもの