フーリエ変換のデータの補間について

質問

Excelでフーリエ変換をする際、データ数は２のべき乗でなければならないと学校で習いました。
データ数が２のべき乗でない場合は２のべき乗になるようにデータを増やす、または減らす必要があるらしいのですが、どうやったらデータ数を増やす、または減らせるのでしょうか？？

例えば、255×255のデータを256×256にするにはどうやったらいいのでしょうか。（MRI画像などは256×256だそうです）
教えてください。

回答No.5

stomachman · ベストアンサー

> MRI画像再構成についての課題です。

　なるほど、それなら話は通じます。実務じゃないんですね。

　周波数空間(いわゆるK空間)の中に、直流成分(空間周波数(0,0))を中心として、等間隔の格子点(n,m) (-N＜n＜N, -N＜m＜N)における複素数のサンプル値s(m,n) = p(m,n)+iq(m,n)が与えられている(p,qは実数)。これを複素逆フーリエ変換して画像を作るためにFFTを使う。

> コンピュータのサンプリングは奇数個でなければならない

　これはおかしい。奇数個でなければならないなんてことはないし、N-1まででやめる理由もない。例えば(-N≦n＜N, -N≦m＜N)の範囲を測ったって良いのですから。しかし、文句を言っても始まらない。ま、この課題に限っては、話の前提条件が(-N＜n＜N, -N＜m＜N)ということなんだな、と考えるしかありません。

> FFTをする場合は２のべき乗個なければならない。

　既に何度も出ている通り、これは数学的には誤りですが、工学的には「さもありなん」というところです。

　ご参考までに、フーリエ変換の定義を少し調べれば分かる、重要な性質を幾つか挙げておきます。

● FFTと逆FFT(IFFT)は、定数倍の違いを無視すれば全く同じ計算。（FFTを2回繰り返すとどうなるか、式を調べてみると面白いでしょう。）

● 2次元FFTは、格子点のそれぞれの行について1次元FFTして、その結果を今度はそれぞれの列について1次元FFTすれば良い。

● 複素IFFTをやるから、結果も複素画像になります。その実部をf(x,y)、虚部をg(x,y)と書くことにすると、
f(x,y) = (1/2)IFFT[(p[m,n]+p[-m,-n])+ i(q[m,n]-q[-m,-n])]
g(x,y) = (1/2)IFFT[(p[m,n]-p[-m,-n])+ i(q[m,n]+q[-m,-n])]
という関係が成り立つ。（従って、g=0となることが分かっている撮影シーケンスの場合（大抵そうです）には、0≦mの部分だけしか測定しないことによって撮影時間を半減することも可能。）

● １次元IFFT（周期M）に喰わせる複素データの列をa[0], a[1], ..., a[M-1]とすると、a[0]は直流成分、a[1]とa[M-1]が周期Mのsin波形、a[2]とa[M-2]が周期M/2のsin波形…の振幅と位相を決める係数です。周期の区切りが0～M-1になってる訳です。
　一方、K空間の方は区切りが-N～Nになっている。FFT(IFFT)はデータが周期Mを持つことを前提にしていますから、
s(m,n)=s(m+jM,n+kM) (j,kは任意の整数）
が成り立つものと考えます。従って、
s(M-m,M-n)= s(-m,-n)
s(M-m,n)= s(-m,n)
s(m,M-n)= s(m,-n)
です。言い換えると、添え字はmodulo Mで0～M-1に書き換えてやる、ということです。

　じゃ本題。N=127、M=256の場合を考えると、s(128,n)(n=-N+1～N-1)やs(m,128)(m=-N+1～N-1)のデータが存在しないからナントカする必要がある。さて、これらは一体何かというと、最も高い周波数の空間周波数成分、つまり画像上で言えば周期2のsin波形（いろんな方向の波がありますが）の振幅と位相を表している。でもその成分については（そして、それより高い空間周波数の成分についても）測定しなかった。

　ならば、そんな成分は画像にも存在すべきでないだろう、というわけで、ここんとこには0を詰める。と考えるのが普通でしょう。
　ところがですね、被写体がシャープでコントラストの大きい輪郭を持っていたりすると、s(127,n)やs(m,127)やs(129,n)やs(m,129)の絶対値が無視できない大きさになる[ANo.4の「両端点」の話です]。すると、0を詰めたことでK空間のs(128,n)やs(m,128)の箇所に鋭い段差が生じます。このような段差は画像上では周期2のsin波形のさざ波(ringing）として現れる。特に、被写体の輪郭に沿って波紋のようなartifactが生じます。
　いや、実は話は逆で、そういう被写体を撮影する以上、正確な画像を作りたかったらもっと高い周波数までMRIで測定すべきなのに、測ってない訳です。

　で、そのringing現象をどう考るか。
「高い空間周波数成分を持つ被写体を撮影したのに、高い周波数については測ってないんだからringingはしょうがない。単に0を詰めとけ。」という「正論」の他に、
「いや、被写体に実在していないringingが出るのはまずい。元来、高い空間周波数成分を持つ被写体を撮影したのに、高い周波数については測ってないんだから、ぼけた画像しか得られないのが当然だ。だからs(127,n)がほぼ0になるようにK空間のサンプルに平滑化フィルタ（|n|と|m|が小さいときはほぼ1で、|n|や|m|がNに近い時には0になるような滑らかな実数値関数h(m,n)）をかけ算しろ。」も実に尤もな話だし、
「確かにringingは困るが、せっかく127の所まで測ったんだからぼかさないでなんとか活かそうよ。少々嘘をつくことになるが、たとえばs(128,n)の値はs(127,n)とs(129,n)の平均値で代用すれば段差は生じない。」という考え方もある。
　ここんとこの判断は、「何を撮影して、結果の画像をどう使うつもりなのか」にも依るでしょう。そういうことまで考えて、いろんな処方について利害得失を検討するのが、ですから、この課題の真髄かと思います。

お礼

クローズが遅れて大変申し訳ありませんでした。
最も詳しく説明してくださった方にベストアンサーをつけさせていただきます。

回答No.4

stomachman

Excelでやると仰っているから、課題か何かでちょっと試している程度のことなのでしょう。でも、もし本気でFFTを使うのであれば、どんな格好のデータを扱っていて[特に両端点はどうなっていて]、FFTをやったあと何をする積もりなのか、に大いに依ります。（なので、もし実務でお困りになってのご質問なら、もう少し詳細を補足して下さい。）

　でもま、大抵の応用では、255点のデータがあったら（256点に増やすどころか）最低でも512点にします。不足したところには0を詰めるの(zero-padding)が普通ですが、そうしない方が良い場合だって多々あります。（たとえば画像において、縁の方の画素が0から遠く離れた値であるなんて時には、安易に0を詰めちゃ駄目な応用が多いです。）

　FFTではデータが周期的であることが前提です。しかし、本来周期なんかない255点のデータのフーリエ変換（従って厳密な意味ではFFTは利用できない）を考えている場合でも、実用上は、うんと長い周期（たとえば512とか1024）の周期関数だと思って扱って構わない場合が多いんで、その時には0(など)を詰めてやることによって長い周期のFFTを利用できるようにする。
　逆に言えば、本来周期なんかない255点のデータを256点のFFTで扱うと、「次の周期」からの干渉を受けて、本来の意味でのフーリエ変換とはだいぶ違う結果になる。なので、どのぐらい以上の周期なら、応用目的に於いて「周期関数だと思って扱って構わない」かを考慮して、何点のFFTを使うかを決めるんです。
　なお、周期が幾らであっもFFTは可能ですけれども、周期が2の冪乗のFFTはプログラムがとびきり簡単に書ける。なので、「うんと長い周期」でありさえすればいいという応用ならば、だったら2の冪乗でいいじゃん、ということで2の冪乗のFFTが使われるんです。

補足

回答ありがとうございます。
確かに、MRI画像再構成についての課題です。

課題内容は
「コンピュータのサンプリングは奇数個でなければならないが、FFTをする場合は２のべき乗個なければならない。どのようにして足りない分を補うのか、どこかの行列からデータを持ってくるのであればどういう理由でその行列を、どこに補うのか（中心部なのか、256列目なのか、それ以外なのか）答えなさい。」
というものです。

自分としてはコンピュータのサンプリングが奇数個という部分からすでに少し意味がわからないです。

回答No.3

sinisorsa

まず既に説明がありましたように、離散フーリエ変換のデータ数は
２のべき乗でなくても、本来定義されます。高速アルゴリズムでは、
２のべき乗の場合に一番効率がよいのですが、Ｅｘｃｅｌでなければ、
もっと自由度があります。
ご質問のご趣旨はＥｘｃｅｌでということなので、これ以上は申し上げません。
さて、本題に戻りますが、
時間信号の場合を例として記します。
フーリエ変換したい信号（データ）は、時間的に有限でしょうか。
（１）有限ならば、フーリエ変換のデータ数というのは、基本帯域
　　　（サンプリング周波数の２分の１まで）の分割数ですから、
　　　足りないデータの分は０を詰めておけばいいのです。
（２）周期的信号の場合には、本来、データ数はその周期内の
　　　サンプル数に一致させる場合がありますが、多くの場合
　　　３周期分程度の有限持続時間信号として、適当な窓関数
　　　（ハミング窓など）を掛けてから、足りないデータとしては、
　　　０を詰めて行います。ただし、線スペクトルとはなりません
　　　ので、分析後に解釈が必要になります。

ご質問は２次元の静止画像のようですね。
この場合、画像は周期的とは考えませんね。
有限範囲に制限されたデータと考えて、
単純に足りないところに０を詰めればよいと
思います。ただし、空間周波数の刻みが少しずれます。
しかし、もともと空間周波数は連続であるので、
問題ないと思います。
試してみてください。

回答No.2

arrysthmia

工学カテゴリーの質問ですね。
数学的には、データ数が 2 の巾でなくても
離散フーリエ変換は行えるし、
データ数を増やしたければ、
0 を与えておけば十分です。

2 の巾乗 という条件は、フーリエ変換を行うのに
FFT という計算手順を使うためのものかと思います。

工学応用に配慮して、
継ぎ目で高周波成分が大きくならないように、
データ数学を増やすには、
滑らかなデータをでっち上げてやる必要があります。

欲しいデータ数周期で周期拡張したデータを
多項式補間する… とか、どうでしょうね。

回答No.1

spring135

フーリエ変換は時間に対しても空間に対しても行われますが、私の経験は時間に対するフーリエ変換だけです。その経験で言うと、255という数字はある時間幅（0～T)を254個の等間隔ΔTに切っていることになります。この場合、同じ時間間隔ΔTで255番目の時間を付け加え256番目のデータを付け加えます。このとき256番目のデ－タの取り方が問題ですが、255番目のデータに比べ、差が大きいと、急に折れ曲がったような分布になり、周波数分布に影響を与えるため、好ましくありません。255番目のデータが0に近ければ256番目のデータを0としてもほとんど影響はありません。明らかに0でない場合は255番目のデータと同じにするか、254番目と255番目のデータの増減を見て、その先の点として滑らかになるような値を選びます。