ママのスキンケアのお悩みにおすすめアイテム

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

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

このQ&Aに関連する最新のQ&A

A 回答 (5件)

> 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をする場合は2のべき乗個なければならない。

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

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

● 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の部分だけしか測定しないことによって撮影時間を半減することも可能。)

● 1次元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)の平均値で代用すれば段差は生じない。」という考え方もある。
 ここんとこの判断は、「何を撮影して、結果の画像をどう使うつもりなのか」にも依るでしょう。そういうことまで考えて、いろんな処方について利害得失を検討するのが、ですから、この課題の真髄かと思います。
    • good
    • 0
この回答へのお礼

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

お礼日時:2010/11/10 16:04

 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をする場合は2のべき乗個なければならない。どのようにして足りない分を補うのか、どこかの行列からデータを持ってくるのであればどういう理由でその行列を、どこに補うのか(中心部なのか、256列目なのか、それ以外なのか)答えなさい。」
というものです。

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

補足日時:2009/06/03 17:54
    • good
    • 0

まず既に説明がありましたように、離散フーリエ変換のデータ数は


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

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

工学カテゴリーの質問ですね。


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

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

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

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

フーリエ変換は時間に対しても空間に対しても行われますが、私の経験は時間に対するフーリエ変換だけです。

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

このQ&Aに関連する人気のQ&A

お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!

このQ&Aを見た人はこんなQ&Aも見ています

このQ&Aを見た人が検索しているワード

このQ&Aと関連する良く見られている質問

Q高速フーリエ変換でデータ数が2のべき乗でない時

こんにちは。現在、フーリエ変換について勉強しているのですが、ちょっとわからないことがあったので質問させていただきました。

質問内容は高速フーリエ変換についてで、cooley&tukeyのアルゴリズムを利用すると、データが2の冪乗個のときは計算量をО(NlogN)に減らせる事ができるというものでした。

しかしデータが2の冪乗個でないとき。例えばN=5000くらいのときはデータを切り取って無理やりN=4096(=2^12)みたいな感じにすれば良いんですよね?
やっぱりその時って、N=5000で通常の離散フーリエ変換したときと周波数値に誤差が出ると思うのですが、それはどうやったら計算できるのでしょうか。。。

どなたかご教授していただければ幸いです。

Aベストアンサー

離散フーリエ変換は、信号が周期的であることを前提としています。
離散フーリエ変換でのデータ数Nは、離散時間信号の周期に当たります。変換の結果は線スペクトルとなります。
N=5000がその信号の1周期なのでしょうか。
もしそうならば、4096にすれば、誤差が大きくなるでしょう。
N=5000で変換すべきです。この場合にも高速アルゴリズムが
存在します。#1の方のとおりです。
FORTRANの時代には、パッケージがありました。
NはN=2^m*3^n*5^k*7^Lだったと思います。

もうひとつの考え方は、有限持続時間信号のフーリエ変換としての
適用です。これは、連続スペクトルとなります。データ数Nは
スペクトルの分解能に関係します。サンプリング周波数をNで割った
ものが周波数分解能となります。
実際のデータよりも2倍程度のNを使うことが多いと思います。
データ数が5000ならば、Nは8192とし足りないデータには、
0を詰めます。これならば、2のべき乗のNを選べます。
この場合、逆変換は周期的な拡張が行われることに注意が必要です。

離散フーリエ変換は、信号が周期的であることを前提としています。
離散フーリエ変換でのデータ数Nは、離散時間信号の周期に当たります。変換の結果は線スペクトルとなります。
N=5000がその信号の1周期なのでしょうか。
もしそうならば、4096にすれば、誤差が大きくなるでしょう。
N=5000で変換すべきです。この場合にも高速アルゴリズムが
存在します。#1の方のとおりです。
FORTRANの時代には、パッケージがありました。
NはN=2^m*3^n*5^k*7^Lだったと思...続きを読む

QExcelで4096点以上のFFTの方法

Excelでは4096点までのFFTしかできませんが、4096点以上のFFTをかける方法(VB?)をご存知の方いらっしゃいましたらお教えください。

Aベストアンサー

下記などご参考まで。プログラムをそのまま公開されています。
http://tsuyu.cocolog-nifty.com/blog/2007/03/publi.html


#プログラムの始めのn=の行とDimの行の数字を適切に修正して使います。

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

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

Aベストアンサー

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

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

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

「時系列の中央部分がゼロに近い」の「中央部分」とは、「IFFTをしたときに得られる、t=0からt=周期までサンプル値を並べた配列の中央付近」ってことですね。周波数空間で応答関数を表現したものが滑らかであればそうなる。応答関数もそれをIFFTしたものも、持っている情報量は同じなので、こりゃしょうがないです。補間をしないで(半分の点数で)IFFTして得られるのは時系列のエリアス(つまり、時系列を周期の整数倍だけ時間軸上でずらしたものを全部の整数について合計したものの一周期分)であって、それ以上...続きを読む

Q周波数解析をエクセルツール『フーリエ解析』で実施したいのですが・・

周波数解析をエクセルツール『フーリエ解析』で実施したいのですが・・
時刻歴を持った波形をエクセルのフーリエ解析にかけたところ以下のような結果が出ました。

時刻(s)波形(風速)(フーリエ結果)
0 0 4964.547892
0.01 0.016148 162.934386092482+757.485796541738i
0.02 0.024223 -176.671853947744+679.108499109482i
0.03 0.032297 -38.1198577747876+304.999881074942i
0.04 0.040371 -33.7184553866481+283.069540754i
0.05 0.047921 67.2878230094194+268.189434427773i
0.06 0.055471 -5.3340068659851+181.755877831686i
0.07 0.063021 43.7700366445313+188.325607658826i
0.08 0.070571 95.1628904739282+115.524911946043i
0.09 0.078121 79.1610244918527+286.809861194846i
0.1 0.084849 -20.9326884388047+207.486443103952i
0.11 0.091578 -10.0874722910491+176.517532576085i

この場合の結果をパワースペクトル密度にしたいのですが,
一番上が定常分で,事項以降が ωt,2ωt,・・・nωtということに
なって,それぞれの実部^2+虚部^2の平方がスペクトルになると思いますが
これを横軸周波数のおなじみのグラフにするには横軸の周波数はどのよう
に考えればいいのでしょうか。
(たとえば2ωtの項の場合の周波数はいくらになる?)
また,風速のパワースペクトルの単位はどのようになるでしょう。
基本的な質問だと思います。とても恥ずかしいのですが,ご教授いただけますか。
もしかすると根本的な間違いをしている気もします・・

周波数解析をエクセルツール『フーリエ解析』で実施したいのですが・・
時刻歴を持った波形をエクセルのフーリエ解析にかけたところ以下のような結果が出ました。

時刻(s)波形(風速)(フーリエ結果)
0 0 4964.547892
0.01 0.016148 162.934386092482+757.485796541738i
0.02 0.024223 -176.671853947744+679.108499109482i
0.03 0.032297 -38.1198577747876+304.999881074942i
0.04 0.040371 -33.7184553866481+283.069540754i
0.05 0.047921 67.2878230094194+268.189434427773i
0.06 0.055471 ...続きを読む

Aベストアンサー

> これを横軸周波数のおなじみのグラフにするには横軸の周波数はどのように考えればいいのでしょうか。

例えば、256点のデータをエクセルで「フーリエ解析」したとします。結果も256点(#0~#255とします)出てきます。ただし、パワースペクトル(複素数の実部^2+虚部^2)にすると、#1~#127と#255~#129は対称形になっていますから、実際には129点の結果になります。
さて周波数ですが、この#1は

 元のデータ256点の中に、ちょうど1周期現れる周波数

に相当します。サンプリングの定理から、最大周波数は、256点の中に128周期の波で、周期は2点毎の波になります(#128)。

具体的に、例えばデータ間隔が1秒で全256点の場合、#1の周波数は1/256Hz、#2の周波数は1/128Hzになります。

Qフーリエ変換について教えてください

フーリエ変換をすると横軸が時間から周波数になるのはわかったのですが、縦軸が何になるのかわかりません。

一般的に縦軸はなにになるのでしょうか?

また横軸が時間で、縦軸が距離をフーリエ変換したら縦軸は何になるのでしょうか?

よろしくお願いします。

Aベストアンサー

時間関数をフーリエ変換すると結果は、その時間関数の周波数成分が
得られます。スペクトルとも言います。従って、縦軸は、周波数成分です。一般に複素数です。
大きさと偏角による表現もできます。
大きさの方は振幅特性、位相角の方は位相特性と呼ばれます。
画像のように空間座標の上の関数の場合には、フーリエ変換すると
空間周波数成分が得られます。横軸は、空間周波数(2次元)となります。
対象とする関数により結果はそれぞれ意味が異なります。
「一般に何になる」とは言えません。

>横軸が時間で縦軸が距離の場合・・・
フーリエ変換の結果は、距離を表す時間関数の周波数成分です。

フーリエ変換の対象の関数は別に時間関数でなければならないということは
ありません。従って、フーリエ変換の結果は適用する人が解釈(定義)すれば
よいと思います。
たとえば、
時間関数をフーリエ変換し、その結果の絶対値の対数のフーリエ変換を
することもあります。これの結果には、発明者らがケプストラムという名前
をつけています。Cepstrum は Spectrum から作った造語です。

時間関数をフーリエ変換すると結果は、その時間関数の周波数成分が
得られます。スペクトルとも言います。従って、縦軸は、周波数成分です。一般に複素数です。
大きさと偏角による表現もできます。
大きさの方は振幅特性、位相角の方は位相特性と呼ばれます。
画像のように空間座標の上の関数の場合には、フーリエ変換すると
空間周波数成分が得られます。横軸は、空間周波数(2次元)となります。
対象とする関数により結果はそれぞれ意味が異なります。
「一般に何になる」とは言えません。

>横軸が...続きを読む

QEXCELにてローパスフィルタを作成する

実験の測定データをEXCELでデータ整理しようと考えております。データ整理のためローパスフィルタをかけたいのですが、具体的にどういった式、もしくはEXCELの機能を使用したらいいのでしょうか?デジタルフィルタが良く分からないのでよろしくお願いします。
ちなみにローパスフィルタは1000Hzをかけたいです。

Aベストアンサー

時系列データの処理ならば

OutputData(n+1) = OutputData(n) + (InputData(n+1) - OutputData(n)) * dt / T

dt:データのサンプリング間隔
T:フィルタの時定数 1/2πf
f:カットオフ周波数
n,n+1:それぞれn個目,n+1個目のデータをしめす。

でいけると思いますが、一次のパッシブなんで効果が薄いかも。(普通はベッセルかけるんでしょうけど、そこまではわからない)

Q窓関数(方形窓)について

EXCEL上でFFTをする際に窓関数を使用したいのですが、方形窓について教えてください。これは取り込んだデータにすべて1をかければよいのでしょうか?あと、方形窓を使用すると周波数分解能があがる?というようなことが調べたらのっていたのですが、どういうことなのでしょうか?どなたかお答え下さい。

Aベストアンサー

というか方形窓は FFT の基本ですよね。

FFT の考え方は、2べき数のサンプリングされた時系列データがあって、それが無限に繰り返されているとしたときの、周波数スペクトル(とその位相)を計算するものです。

(おっと、本来は逆ですね。以下、蛇足です。無限に繰り返される時系列データを、方形の窓関数で切り出して、スペクトルを計算する、というのが離散フーリエ変換です。ついでながら FFT とはちょうどサンプリングデータ数が2べき数の時に素早く計算できるアルゴリズムのことで、基本的な性質は離散フーリエ変換と同じです。)

ところがこれをそのまま実際のデータに応用すると、サンプリングの最後のデータの次のデータは(時系列は無限に繰り返されるので)、最初のデータになります。そのときにちょうどそのような感じでデータがたまたまサンプリングできれば良いのですが、実際そうはうまくいきません。で、どうなるかというと、最後のデータと最初のデータが急に大きく変化するような具合になって、その影響が FFT の結果に現れてます。つまり本来ない周波数成分が、最初と終わりの不連続(?)なデータによって現れたりします。

この問題を回避するのが方形以外の様々な「窓関数」です。窓関数を見ればお分かり頂けるとおもいますが、いずれも最後と最初のデータを小さくするようになっており、上記の影響を少なくしようとしているのが分かると思います。しかしながら、窓関数の影響も勿論あって、その多くの場合、得られるスペクトルがピークが低く、幅が広がるように「ぼけた」感じになってしまいます。この意味で、質問者さんがおっしゃるように「周波数分解能が下がる」わけです。

というか方形窓は FFT の基本ですよね。

FFT の考え方は、2べき数のサンプリングされた時系列データがあって、それが無限に繰り返されているとしたときの、周波数スペクトル(とその位相)を計算するものです。

(おっと、本来は逆ですね。以下、蛇足です。無限に繰り返される時系列データを、方形の窓関数で切り出して、スペクトルを計算する、というのが離散フーリエ変換です。ついでながら FFT とはちょうどサンプリングデータ数が2べき数の時に素早く計算できるアルゴリズムのことで、基本的な性質は離...続きを読む

Qエクセルを使用してデジタルフィルタのハイパスフィルタ・ローパスフィルタの掛け方を教えてください。

タイトルどおりなのですが、エクセルを使用してデジタルフィルタのハイパスフィルタ・ローパスフィルタをかけたいです。

ご存知の方、よろしくお願いします。

Aベストアンサー

わたしはわかりませんが、こんなページがありました。

エクセルを使用してデジタルフィルタのハイパスフィルタ・ローパスフィルタの掛け方を教えてください。
出来ましたらそのまま使える式をお願いいたします。

http://q.hatena.ne.jp/1204546061

Qフーリエスペクトルの振幅について

ある時間関数を離散フーリエ変換して得られるフーリエスペクトルの
振幅値について教えて下さい。

今想定している離散フーリエ変換の式は一般的なもので
Σ(k=0~N-1) f(k)exp(-2πkni/N)
を考えています。
また、離散フーリエ変換して得られるスペクトルは
√(Re^2+Im^2)
で計算します。

離散フーリエ変換を適用する関数を、
振幅1の直流、及び振幅1で周波数5[Hz]の正弦波とします。
(この2つの信号は別々の信号で合成されていません。)
サンプリング周波数を20[Hz]とした場合、
サンプリングして得られるデータ列はそれぞれ、
直流: 「1, 1, 1, 1」
正弦波: 「0, 1, 0, -1」
となると想定されます。
(正弦波をサンプリングする場合は位相が関わってきますが、
今回は気にしないで下さい。)

このデータ列に対して上記の離散フーリエ変換を適用した場合、
得られるフーリエスペクトルの振幅値はそれぞれ、
直流: 「4」(直流のフーリエスペクトルの振幅値値)
正弦波: 「2」(5[Hz]のフーリエスペクトルの振幅値)
となります。
(データ点数は上の通り4点)

ここで質問なのですが、
離散フーリエ変換して得られるスペクトルの振幅値から
元の関数の振幅値を求める場合、
フーリエスペクトルをサンプリングの総データ点数で割ることは
数学的に納得できます。
しかしこの例の場合、フーリエスペクトルを総データ点数で割ると、
直流: 「4 -> 1」
正弦波: 「2 -> 0.5」
となってしまい、直流は正しいのですが、
正弦波の元の振幅値を正確に求めることは出来ません。
フーリエスペクトルの振幅値から正弦波の振幅値を正しく求めるには、
「フーリエスペクトルの振幅値*2/データ点数」
としてやらなければいけません。

上記のことに関して、
なぜこのようになるのかを(2をかける理由を)教えて頂けないでしょうか。

ある時間関数を離散フーリエ変換して得られるフーリエスペクトルの
振幅値について教えて下さい。

今想定している離散フーリエ変換の式は一般的なもので
Σ(k=0~N-1) f(k)exp(-2πkni/N)
を考えています。
また、離散フーリエ変換して得られるスペクトルは
√(Re^2+Im^2)
で計算します。

離散フーリエ変換を適用する関数を、
振幅1の直流、及び振幅1で周波数5[Hz]の正弦波とします。
(この2つの信号は別々の信号で合成されていません。)
サンプリング周波数を20[Hz]とした場合、
サンプリングし...続きを読む

Aベストアンサー

離散フーリエ変換というのは、実質離散フーリエ級数展開なので、
フーリエ級数展開を考えます。

f(t) = a0/2 + Σ[n=1→∞] { an cos (nωt) + bn sin (nωt) }

f(t)の周期をTとして、ω=2π/Tです。

直流成分の振幅といっているのはこの第1項a0/2のことで、
サイン成分の振幅はbnのことです。

問題文の離散フーリエ変換の式

>Σ(k=0~N-1) f(k)exp(-2πkni/N)

は複素フーリエ級数展開なのでオイラーの公式

cos nωt = [e^{inωt}+e^{-inωt}]/2
sin nωt = [e^{inωt}-e^{-inωt}]/2i = = -i [e^{inωt}-e^{-inωt}]/2

を使って書き直すと、

f(t) = a0/2 + Σ[n=1→∞] { ([an -i bn]/2) e^{nωt} + ([an +i bn]/2) e^{-inωt} }

an = a(-n), bn = -b(-n)の関係があるので、

cn = ([an -i bn]/2), c(-n) = ([an +i bn]/2)

と置くことができ、a0/2をc0と定義し直せば、

f(t) = Σ[n=-∞→∞] cn e^{nωt}

したがって、複素フーリエ係数が求めているのはcn = (an-ibn)/2で、その実数部はan/2、虚数部は-bn/2です。

こうなる理由は、サイン、コサインのときは正の整数だったnを複素数で取り扱うときにマイナス側に拡張したことで、同じ係数が+側と-側にわかれたためです。

離散的な場合は和が-N/2~N/2の範囲の有限項で打ち切られ、
-N/2~0の範囲が一周期ずらされてN/2~Nになっています。

離散フーリエ変換というのは、実質離散フーリエ級数展開なので、
フーリエ級数展開を考えます。

f(t) = a0/2 + Σ[n=1→∞] { an cos (nωt) + bn sin (nωt) }

f(t)の周期をTとして、ω=2π/Tです。

直流成分の振幅といっているのはこの第1項a0/2のことで、
サイン成分の振幅はbnのことです。

問題文の離散フーリエ変換の式

>Σ(k=0~N-1) f(k)exp(-2πkni/N)

は複素フーリエ級数展開なのでオイラーの公式

cos nωt = [e^{inωt}+e^{-inωt}]/2
sin nωt = [e^{inωt}-e^{-inωt}]/2i = = -i [e^{inωt}...続きを読む

Q不規則なデータのfft処理

データサンプリングが不規則なのですが、ExcelでMacroを使って中間値を補間してフーリエ変換FFT処理する方法をお教え頂けないでしょうか?
またできればExcelはFFTできる最大が4096ですがMacro処理でそれ以上の数でできるようになりませんか?フリーソフトでも構いません。切り貼りは致します。
両方の質問共にでなくても、片側だけでも良いので、お手数おかけしますが、よろしくお願いいたします。

Aベストアンサー

マクロが組めるならば、無理にFFTをせずに、DFT(離散フーリエ変換)で良いのかと思います。
また、かなり誤解がある(一般にそうだろうと思いますが)ようなので、
少し難しくなりますが、頑張って読んでください。

FFTと言う名前は有名ですが、DFTと演算結果は同じです。
DFTの演算量を減らして高速化したものがFFTです。

DFTは非常に簡単なループ処理で演算できます。
FFTは結構面倒なので、一からプログラムすると大変です。

私はFFTのサンプルを手に入れるのが面倒なので、
必要に迫られたときネットで検索してDFTの数式を見つけて毎回プログラムしてます。

【これを知っておくと、C言語だろうが、Javaだろうが、エクセルだろうが、何でも移植できる】

プログラムも総ライン数で10行くらいでしょうか。
デジタルデータの周波数変換に関してFFTを頼る必要は全く無いですよ。

【DFTで検索。Σ式がすぐ出てきます。つまり数列の和ですね。表計算で簡単に出力できるでしょう】

FFTは今から相当昔に開発された手法ですが、コンピュータが遅い時代に重宝しました。
しかし、当時のマイコンでもDFTで充分なスピードでしたよ。
今のパソコンのスピードは、当時の1000倍くらいです。

定義をおさらいします。

(1)フーリエ変換
 無限数列を前提に処理するもの。周波数変換。
(2)DFT
 サンプリングされたデータを擬似的に周波数変換するもの。厳密に周波数変換ではない。
(3)FFT
 DFTと演算結果が同じになるように開発された高速演算アルゴリズム。
 アルゴリズム名であり、変換としての用語ではない。直交変換としてはDFTが正しい説明。
 DFTの演算結果を見ると、周期的に、実数と虚数の絶対値が同じものが現れたりする事に着眼。
 最初のほうだけDFTを行い、途中から虚実を入れ替えたり、符号を変えたりして、
 残りの演算結果を生成してしまう手法。
 
FFTが有難いとされていたのは、(3)の説明で分かると思いますが、浮動小数点演算よりも、
符号の入れ替えや、実数部と虚数部の入れ替えの方が早いとされていた時代の話です。

3Dポリゴンのゲームなどは、浮動小数点を使った行列演算を物凄いスピードでやっております。
と言うことは、DFTとFFTの違いにこだわる必要もなくなってきたわけです。

その上で、DFTは普通の人でもプログラムできます。正に初心者向け。

次にサンプリング周期が不規則とありますが、これはディジタル信号処理の感覚ではNGです。

せめてサンプリングレートを落として、一定周期でサンプリングされていることを保障しましょう。
サンプリングしたタイミングがズレると、歪みの様なノイズが重畳されます。(変調がおきる)

さらに、通常のFFTの使い方は、周波数積でのフィルタリングです。
この用途の一つとしてデータ補完があります。
インターポレーションというテクニックです。
これにより粗いサンプリングをしたデータの中間値を補完する事ができます。

これは直線補完や二次補完よりも滑らかな補完になり、自然に見える効果があります。

アニメやイラスト画像には向いていませんので、スプラインとか色んな補完が使われます。
しかし、音声や自然画像(写真)の補完では周波数をつかった補完が効果的で重宝されます。
JPEGやMPEGなどで使われるDCTもこれと似た考え方をしています。
つまり、一度圧縮でデータを削っておいて、補完によって復元するわけです。


そのため、

【FFTする前に別のアルゴリズムで中間値を補完してしまったら、意味が無いです】

中間値を直線的に案分して補完したり、投影法により補完していたのは相当昔(40年以上前)
の話であり、今はDFTを使って補完するほう普通です。(特性が良い)

また、
DFTでは少ないサンプルデータでも大雑把に周波数のカーブを見せてくれるので便利なのです。

例1)サンプルデータは5個しかないが、周波数データは1024サンプル欲しい。

【1024の配列を用意し、先頭5個にサンプルデータを代入。残りに0を入れてDFTをする】

 これでOKです。

例2)データサンプルを補完を使って増やしたい(インターポレーション)
 元データ:A、B、C、D、 の4サンプルを例に、データを2倍に増やす方法を説明します。
 手順1) A、0、B、0、C、0、D、0 と言う0(ゼロ)を間に挿入したデータを作る。
 手順2)DFTをかける
 手順3)DFTデータを見ると、周期的に同じデータ出てくるので2周期目を虚数部含めて
     0(ゼロ)で消す
     a,b,a,b,b,a,b,a 見たいな感じに成るので、a,b,0,0,0,0,b,a に作り直す
 手順4)IDFT(DFTの逆変換)をする。これも式が違うだけでDFTとほぼ同じ演算方法です。
 手順5)IDFT結果を見て、元のデータと似たような変化を持つ数列になっていることを確認する。

例3)歯抜けデータを補完する
 歯抜けになった部分を0として例2と同様の処理をする。
 この場合サンプリングレートは半分になります。また完全な復元ではありません。
 大雑把に復元できた気分になります。

以上ご参考になれば。

マクロが組めるならば、無理にFFTをせずに、DFT(離散フーリエ変換)で良いのかと思います。
また、かなり誤解がある(一般にそうだろうと思いますが)ようなので、
少し難しくなりますが、頑張って読んでください。

FFTと言う名前は有名ですが、DFTと演算結果は同じです。
DFTの演算量を減らして高速化したものがFFTです。

DFTは非常に簡単なループ処理で演算できます。
FFTは結構面倒なので、一からプログラムすると大変です。

私はFFTのサンプルを手に入れるのが面倒なので、
必要に迫られたときネットで検...続きを読む


人気Q&Aランキング