FFTにおけるゼロ追加、補間や分解能について

質問

FFT解析で、時間軸データの後にゼロをつけると、周波数分解能が上がることは数値計算で確認しましたが、次の場合のやり方がわかりませんので教えていただけないでしょうか。
1.　周波数軸データの後にゼロをつけて最高周波数をあげ、逆変換して刻み時間を細かくしたいのですが、ゼロの配置は共役関係やナイキスト周波数を考えるとどういう風になるのでしょうか。
2.　周波数軸データの値を内挿すると周波数分解能が上がり、逆変換すると結果として時間軸の継続時間が長くなるように思うのですが、それでよいのでしょうか。またその場合、良い内挿の方法はあるのでしょうか。
3.　以上に関するURLや書籍を紹介していただけないでしょうか。

回答No.4

stomachman · ベストアンサー

「時系列の中央部分がゼロに近い」の「中央部分」とは、「IFFTをしたときに得られる、t=0からt=周期までサンプル値を並べた配列の中央付近」ってことですね。周波数空間で応答関数を表現したものが滑らかであればそうなる。応答関数もそれをIFFTしたものも、持っている情報量は同じなので、こりゃしょうがないです。補間をしないで（半分の点数で）IFFTして得られるのは時系列のエリアス（つまり、時系列を周期の整数倍だけ時間軸上でずらしたものを全部の整数について合計したものの一周期分）であって、それ以上のことは情報がないから分からない。
　もしも何か別の拘束条件があるのなら、一応、手がない訳ではありません。すなわち、拘束条件によって時系列関数としてあり得る関数が絞られるから、その中から解を探すんです。しかしこの手（「超解像」という）は誤差に極めて敏感で不安定であり、だから応答関数の精度が非常に良いと保証できる場合でない限り使えません。（10時間掛けて計算するぐらいだから、そんな保証はタブンできないでしょう。）拘束条件がないのなら、せいぜい例えば「FFTしてサンプルを間引いたときに応答関数に一致するような時系列のうちで、ノルムが最小であるもの」を「一例として」構成してみせるぐらいが関の山です。

　なお、もしもご質問が「時系列のうち、後端ならともかくヨリニヨッテ中央部分が抜けるのが困る（「後端部分はどうせゼロに近い」とアプリオリに分かっていて、だからIFFTでドンナ値が出ようと誤差に過ぎないと考えて利用しない）」という話であるのなら、それは時間軸の原点の位置がまずいだけかも知れません。すなわち、ひょっとすると「応答関数を細かく計算するとギザギザしていて、それを１点おきに間引いてから補間したためにギザギザが失われ、IFFTすると時系列が時間軸上で周期/2だけずれて現れる」という現象なのではないか。
　どういうことかと言いますと、ある滑らかな応答関数をIFFTして得た時系列（中央部分がほとんど0）を考え、その時系列を時間軸上で周期/2だけ平行移動したとします。それをFFTして得られる応答関数はギザギザに振動していてでも１点おきに間引くと元の滑らかな応答関数の格好になる（。で、これを10時間掛けて計算したのだとしましょう）。この、間引いた物（滑らかな応答関数）を普通のヤリカタで補間して得られるものも当然滑らかであり、これをIFFTすると、最初の時系列（中央部分がほとんど0）、つまり時間軸上で周期/2だけずらす前の物に（ほぼ）戻ってしまう。
　もしそういうことが起こっているのだとすると、計算なさった時系列の後端部分はゴミなんかじゃなく、それを前端部分の前にくっつけたものこそが時系列の本来の姿ということになります。

　またしても滑ってるかも知れませんが…

回答No.3

stomachman

ANo.1、もういっちょ滑ってました。

>・ 時間軸データの値を内挿する　⇔ 時間軸データの値を内挿する

同じに決まってる。いや、そうじゃなくて、

・時間軸データの値を内挿する　⇔ 周波数軸データの後にゼロをつける

です。

補足

指のすべりの件、多分そうだろうと解釈していました。
一連のご回答を見て、この分野に関して非常にハイレベルの方とお見受けしております。
実は地震波動解析の際に周波数応答関数を計算し、それを逆変換して時系列（地震波)を計算しています。ただ周波数応答関数の計算に10時間前後かかるため、時間短縮のため計算周波数点を粗くし、それを補間してから逆変換で時系列を求めたいのが目的です。実際に計算してみると、時系列の中央部分がゼロに近い値となり、目的の時系列が得られないようです。
これはご指摘の、
・周波数軸データの値を内挿する　⇔ 時間軸データの後にゼロをつける
の表れと解釈しておりますが、この現象が生じないように、目的の計算結果を得る方法はないものでしょうか。

回答No.2

stomachman

ANo.1でうっかり滑りました。

一番最後の「こいつを逆フーリエ変換したものが」の「こいつ」ってのは、矩形窓を掛ける前の「ΣA(ω+2kΩ)」のことで、逆フーリエ変換の結果はb(t)に他なりません。

補足

前回の補足説明、訂正のつもりが誤って投稿ボタンを押してしまいましたので、注の部分を以下に訂正します。
注：理論的に詰めた訳ではなく、直感的に次のように考えた。
元のデータではナイキスト周波数成分は1個であるが、ゼロ追加後の
データでは共役をとる関係でこの成分は2回表われるので、値を半分にする。

回答No.1

stomachman

FFTとその逆変換は、軸の正負が逆になることと定数倍以外、全く同じ変換です。そして、
・時間軸データの後にゼロをつける　⇔ 周波数軸データの値を内挿する　
・時間軸データの値を内挿する　⇔ 時間軸データの値を内挿する
　内挿をsinc関数の畳み込みによって理想的に行った場合、「⇔」で結んだ一対の操作は同じことです。

　数値的に色々試していらっしゃるようですけど、フーリエ変換と離散フーリエ変換（FFT)の関係をある程度把握すると、「連続関数をサンプリングする」「補間する」などの操作が全てフーリエ変換の言葉で表現できるようになり、見通しが良くなると思います。こんな具合です：

　ある離散データのサンプルひとつf[k]を、δ関数 f[k]δ(t-kΔt)だと考え、それらの総和
　g(t)=Σf[k]δ(t-kΔt)　(Σはk=-∞,…,∞）
をひとつの「連続データ」だと思うことにします。なんでこんなことするかと言いますと、離散データをこういう形に表すことによって、その（FFTではなく）フーリエ変換を考えることができるからです。
　f(t)をフーリエ変換すると、FFTで得られるスペクトルと同じものが周波数軸上に周期的に無限個並んだ形G(ω)が現れます。その周期Ωはサンプリング間隔Δtの逆数に比例します。
　そして、G(ω)の波形から周期ひとつ分だけを取り出して、逆フーリエ変換して得られるのは、元の離散データg(t)にsinc関数を畳み込んで作った連続データです。なぜなら、「sinc関数を畳み込む」のは、「フーリエ変換し、矩形窓をかけ算し、逆フーリエ変換を行う」のと同値。そして、「周期のうちのひとつ分だけを取り出す」操作は「矩形窓をかけ算する」ことに他ならないからで、その矩形窓の幅Ωはサンプリング周波数（ナイキスト周波数の２倍）です。

　次に、離散データに対して一次補間を使って、サンプリング間隔を半分にした離散データを作る場合を考えましょう。
　一次補間とは二等辺三角形をした関数s(t)を畳み込むことに他なりません。この三角形の底辺の長さはサンプリング間隔Δtの丁度2倍で、これは「サンプリング間隔Δtの幅を持つ矩形」二つを畳み込んだ物と同じです。したがって、s(t)のフーリエ変換は「サンプリング周波数Ωのところに最初の零点があるsinc関数」の2乗。これをSとしましょう。
　で、先に述べた「同じものが周波数軸上に周期的に無限個並んだ形」GにSをかけ算し（従って周期関数ではなくなります。これをA=SGとしましょう）、さらに逆フーリエ変換したのが一次補間で作った連続データ（折れ線グラフになりますね）a(t)です。
　次に、a(t)を元のサンプリング間隔Δtの半分の間隔Δt/2で再度サンプリングし（そして無限個のδ関数の総和で表し）たものをb(t)とします。これをまたフーリエ変換すると、「同じものが周波数軸上に周期的に無限個並んだ形」Bが現れますが、その周期はGやAの周期の２倍(2Ω)になっています。これは、さっき見たAを周波数軸に沿って周期2Ωの整数(k)倍だけずらしたものすべてを総和した波形
ΣA(ω+2kΩ)
に他なりません。で、これに幅2Ωの矩形窓をかけ算して周期のひとつだけを取り出してやると、これがAの「エリアス(alias)」です。こいつを逆フーリエ変換したものが、「一次補間でサンプリング間隔を半分にした離散データ（をδ関数で表したもの）」になっています。

補足

非常にご丁寧なご回答、ありがとうございます。今のところご回答の中身を完全に理解しきれていませんが、当方の質問の1、2については分かりました。質問の1について（恥ずかしながら数値計算で）確認したところ、元のFFTによるフーリエ変換のナイキスト周波数の値を半分にし、ナイキスト周波数より高い周波数成分にゼロを追加後、共役をとって逆変換すれば、元の時系列を補間できました。例えばデータ個数nの時系列を補間して2n個の時系列を得る場合をfortran風に書けば
x(t)⇒X(iω)=X(k), t=1～n
X(n/2+1)=X(n/2+1)/2注
X(k)=0, ｋ=n/2+2～2n
X(n-k+2)=conjg(X(k))
X(k)⇒x(t), t=1～2n
注：理論的に詰めた直感的に次のように考えた。
元のデータではナイキスト周波数成分は1個であるが、ゼロ追加後の
データでは共役をとる関係でこの成分は2回表われる、半分にする。