hirax.net::Keywords::「FFT」のブログ

■Wveletで周期解析をしてみる　　

　音声・地震などの1次元信号や、画像等の2次元信号処理の解析というのはなかなか面白そうだ。そこで、周期ムラに対してWaveletをかけて周波数解析をする練習をやってみたい。また、短時間フーリエ変換とWaveletの比較もしてみたい。音声・地震などのデータはまた別にやってみることにして、今回は画像データを扱うことにする。ただし、いきなり2次元も何なので、画像データの周期(つまり1次元的な振動)に注目して、解析を行ってみたい。

　まずは、「周期ムラのある画像」と「周期ムラのない画像」の2種類の画像を作成する。画像はいずれも数式を用いて作成した。X方向に変化する縞模様であり、表.1のような演算式になっている。一応、2次元画像ではあるが、Y方向にはなんの変化もない。2つの数式を見比べてみると判るが、いずれも2項からなり、低周波数のSinと高周波数のSinからなっている。「周波数ムラのある画像」では、その低周波数のSinの中にさらにSinがあるので、周波数がある周期で変化していることになる。一見、「周波数ムラのない画像」の方でも低周波数のSinの内部にさらにSinがあるように見えるが、0が掛けられているので、実際には存在しないのと同じである。

表.1:画像を作成するために使用した演算式

2つのSinからなり、その一方のSinの周期がムラ(一定の周期)をもっているもの。第二項目のSinの内部にさらにSinを入れることにより、周波数ムラを作っている。	2つのSin波からなり、どちらの周期も正確なもの。第二項目のSinの中のSinは0をかけてあるので、何ら影響を及ぼさない。

　そのような数式に基づいて作成した画像を図.1に示す。なお縦軸がX軸であり、横軸がY軸である。図.2(b)では周波数ムラはないが、2つの周波数成分から作成されているため、うねりが生じている。

図.1:作成した原画像

周波数ムラのある画像	周波数ムラのない画像
図.1(a):	図.1(b):

　まずは、Wavelet変換である。図.2がその結果である。縦軸が周波数を示している。縦軸の上方向が高周波を示し、下方向が低周波を示している。また、横軸が原画像のX方向である。白は強度が小さいことを示し、黒は強度が強いことを示している。
　いずれの画像も2つの周波数成分からなることが一目瞭然である。また、図.2(a):「周波数ムラのある画像」の方では低い周波数成分の方が、さらにある周期で周波数が変化していることがわかる。

図.2:Wavelet解析を行ったもの (Daubechiesの6次のCoifletFilterを使用)

周波数ムラのある画像	周波数ムラのない画像
図.2(a):	図.2(b):

　同じWavelet変換でも異なるFilterを用いてみると、結果は異なる。例えば、図.3がその例である。こちらの方が「周波数ムラ」がどのように生じているかを見るにはいいかもしれない。

図.3:Wavelet解析を行ったもの (Daubechiesの6次のLeastAsymmetricFilteｒを使用)

周波数ムラのある画像	周波数ムラのない画像
図.3(a):	図.3(b):

　それでは、Wavelet変換ではなくて、フーリエ変換を用いて周波数解析を行ってみる。先ほどの1次元データの全領域に対してフーリエ変換をかけてみる。その結果が図.4である。ここで、横軸が周波数を示し、右側が高周波数を示し、左側が低周波数を示している。縦軸は強度である。
　このフーリエ変換の場合も、2つの画像が2つの周波数成分からなり、図.4(a):「周波数ムラのある画像」では低周波数成分がぶれているのはわかる。しかし、その周波数ブレがどのようなものであるかまでは、わからない。

図.4:全領域にFFTをかけて、周波数解析を行ったもの

周波数ムラのある画像	周波数ムラのない画像
図.4(a):	図.4(b):

図.5:短時間FFTをかけたもの

周波数ムラのある画像	周波数ムラのない画像
図.5(a):	図.5(b):

　単なる全領域にわたった周波数解析と、位置と周波数が同時にわかる解析の違いは非常に大きい。使いこなすのはなかなか難しそうだが....

■オシロスコープソフトを作る　

PCを2Ch高性能オシロスコープにしたい

　オシロスコープがあると便利だが、家で使うには敷居が高いし、値段も高い。まして、FFTアナライザーがついて周波数解析も行うことができるような機械になると、遊びで買うという値段ではなくなってしまう。そこで、PCを2Ch高性能オシロスコープにするソフトを作ってみたい。

　以前、音階を調べた時に、SoundBlaster互換I/Fを使ったLabViewのサウンド入出力のViを使用してみる。目標はとにかくオシロスコープと同じ使い勝手であること、使うのが簡単であること、そして、周波数解析などが簡単に行うことができること、である。入力として、音声入力を使用しているので、たいていのPCで使うことができるし、音声入力マイクが着いているPC(たいていのノートPCは着いているだろう)なら、マイク(あるいは入力端子)を用意する必要すらない。

　というわけで、下がそのアプリケーションの画面。

作ったアプリケーションの画面 (口笛を吹いて、周波数解析を行ったところ)

　20KHz(ナイキスト周波数で言えば正確には10KHzか)までしか、使うことはできないが、ちょっと使いたい時には便利だ。特に、音声を解析したいならば、必要十分である。上の画面では口笛を吹いて、その音声波形を表示させ、周波数解析を同時に行っている。周波数ピークが表示され、1120Hzであるという表示がされる。

　内蔵マイクを使用すると、ただアプリケーションを走らせれば、音声周波数解析が行える。もしも、比較的高性能なマイクがPCに着いているならば、リアルタイム振動解析すら行うことができる。もちろん、マイク入力端子に何らかの入力をすれば、どんな解析もできるわけだが、何の用意もせずにできるというのは便利である。例えば、うるさいデスクトップPCの近くへこのアプリケーションを走らせたノートPCを近づけると、デスクトップPCがなぜうるさいかを簡単に調べることができる。うるさいのは、ハードディスクの周波数なのか、ファンの周波数なのかすぐにわかる。
　
　今回、作成したアプリケーションはここにおいておく。動作させたら、レンジを調整することを忘れずに。
Ocilo.lzh LZH形式　1,259kB (打ち止めです。あしからず。)
LabViewのアプリケーションライセンス上、ダウンロード数は50回までで、その数近くになったら削除することにする。

■振動・声紋解析用のソフトをつくる　

PCオシロソフトを高機能にしたい

　前に作った2Chオシロ&FFTアナライザーに時間vs周波数グラフの表示機能を付けたい。そうすれば、もしも音声解析に利用するならば、声紋分析もできる。また、振動解析ならば、周波数変化を簡単に調べることができる。こういったものが簡単に作れるのはLabViewの素晴らしい所だ。

　このソフトを動かすと、コンピューターのファンやハードディスクの回転数はとても鋭い周波数ピークを持っていることがよくわかる。また、マイクに向かってしゃべれば、声紋分析も可能である。ウソ発見器などにも応用してみたい。

　下がアプリケーションを動かした様子である。左上が生波形、左下が周波数vs強度、右上が1Ch目の時間vs周波数分布、右下が2Ch目の時間vs周波数分布である。

作ったアプリケーションの画面 (適当にしゃべったところ)

　測定用の準備は整ったので、このアプリケーションを使って色々な音声解析や振動解析などをしてみたい。ところで、このアプリケーションを作成した所で、ほとんど同じようなソフトの広告を見かけた。

http://www.mcor.co.jp/goods/fft/

　上に載っているソフトと同じようなことは今回のソフトを使えばできると思う。また、http://www.mcor.co.jp/goods/fft/にあるソフトの便利そうなところは参考にしたい。

　ここに今回作成したアプリケーションを置いておく。
ocilo2.lzh 1,266KB (打ち止めです。あしからず。)
　LabViewのアプリケーションライセンス上、ダウンロード数が40本を近くなったところで削除する。

■音場の定位を見てみたい　

立体音感を考える　その2

　前回(といっても間に他の話も挟まっているのだが)、

• 立体音感を考える - バーチャルサウンドソフトウェアを作ってみよう- (1999.12.06)

　色々なことを考えるには、その目的にあった測定器が必要である。何か新しいことをしようと思ったら、そのための新しい測定器を作成しなければならない(と思うだけだが)。そして、何より私は計測器なんてほとんど持っていない。だからといって、計測器を買うお金があるわけではない。というわけで、困ってしまうのだ。

　そこで、立体音感を考えるための測定器を作っていくことにした。といっても、すぐにできるとも思えないので、色々実験をしながらボチボチとやってみることにした。勉強がてら、ボチボチやってみるのである。オーディオ関連のことにはかなり疎いので勉強にはちょうど良いだろう。

　資料をいくつか眺めてみたが、特に

• 「立体視の不思議を探る」　井上　弘著　オプトロニクス社

• 音像定位の因子
• 両耳差因子 (音響信号)
• 音の強さ(振幅)の差
• 位相の差
• 周波数スペクトル因子

　そこで、いきなりだが今回作成した解析ソフト「音場くん一号」のアルゴリズムは以下のようになる。

1. PCのサウンド入力から、サンプリング周波数 22.05kHz、Stereo 各チャンネル8bitで取り込みを行う。
2. 取り込んだデータを4096点毎にウィンドウ(Hamming or無し)処理をかける。
3. 高速フーリエ変換(FFT)を行う
4. FFTの結果の実部について、左右のチャンネルの差分を計算する

　次に示すのが、「音場くん(仮名)一号」の動作画面である。「音場くん(仮名)一号」の画面構成は、

• 右側->制御部
• 左側->計測データ表示部

• 音声波形データ(赤=左、緑=右)
• 周波数(横軸)vs左右での音圧の差(縦軸)
• 時間(横軸)vs周波数(縦軸)vs左右での音圧の差(色)

「音場くん(仮名)一号」の画面
「安売り水着を結局買ったアタシの歌」イントロ部

(黒字に赤、緑の色構成は変更の予定)

　計測データ表示部の拡大図を下に示す。

• 音声波形データ(赤=左、緑=右)
• 周波数(横軸)vs左右での音圧の差(縦軸)
• 時間(横軸)vs周波数(縦軸)vs左右での音圧の差(色)

　この表示計の意味を例を挙げて説明したい。例えば、下の画面では左の方に定位している音が鳴ったときの状態を示している。一番上の音声波形データでは緑(右)の波形は小さいのに対して、赤(左)の大きな波形が見えている。
　また、真ん中の「周波数(横軸)vs左右での音圧の差(縦軸)」では横軸100(任意単位)程度の高さの辺りで左チャンネルに位置する音が発生しているのがわかる。
　また、一番下の「時間(横軸)vs周波数(縦軸)vs左右での音圧の差(色)」では時間的に一番最後(横軸で右側)の方の横軸560、縦軸100位の位置に白い(すなわち左チャンネルに定位する)音が発生しているのがわかると思う。

「音場くん(仮名)一号」の画面の拡大図
「安売り水着を結局買ったアタシの歌」イントロ部

　この曲のイントロでは、「ポンッ」という音が高さを変えつつ、左右にパンニング(定位位置を変化させること)する。
　一番下の「時間(横軸)vs周波数(縦軸)vs左右での音圧の差(色)」を示したグラフ中で白・黄色(左に定位)と青・黒(右に定位)する音が時間的にずれながら現れているのが判ると思う。

　このようにして、この「音場くん(仮名)一号」では音の定位状態についての「極めて大雑把な」計測が可能である(保証はしないけど)。「音場くん(仮名)一号」を使った他の例を示してみる。

　下は種ともこの「O・HA・YO」の中から「The Morning Dew」のイントロ部を示したものだ。

• 左(白・黄)チャンネル方向に定位するピアノ
• 右(黒・青)チャンネル方向に定位するガットギター

「The Morning Dew」のイントロ部での
「時間(横軸)vs周波数(縦軸)vs左右での音圧の差(色)」
を示したもの

　これはまるでオルゴールのピンを見ているようだ。あるいは、シーケンサーや昔の自動演奏ピアノのロール譜のようである。対位法などの効果をこれで確認したくなってしまう。

　さて、ここまでの例は楽器も少なく、比較的自然な定位状態であった。しかし、以下に示すような場合には不自然なくらいの「音の壁」状態の場合である。かなり状態が異なる場合だ。

「KI・REI」のラストのラストコーラス部での
「時間(横軸)vs周波数(縦軸)vs左右での音圧の差(色)」
を示したもの

　これは、種ともこの「O・HA・YO」の中から「KI・REI」のラストのラストコーラス部を示したものである。人のコーラスが重なり合っていく部分である。色々な高さの声が重なり合っていく様子がわかるだろう。
　ところが、このグラフをよくみると、同じ音が時間的に持続しているにも関わらず、時間毎に定位位置が左右で入れ替わっているのがわかる。

　これはきっとエフェクターで言うところのコーラスなどをかけたせいだろう(素人判断だけど)。人工的にフィルタ処理をしているためにこのようになるのだろう。こういう結果を見ると、「音場くん(仮名)一号」をプログレ系の音の壁を解析してみたくなる。

　さて今回は、音声の定位状態を解析する「音場くん(仮名)一号」を作成し、いくつかの音楽に対して使ってみた。まだまだ「音場くん(仮名)一号」は作成途中である。これから続く立体音感シリーズとともに「音場くん(仮名)」も成長していく予定である。

　さて、一番先の画面中に"Re"という選択肢があるのがわかると思う。もちろん、これと対になるのは"Im"である。FFTをかけた結果の"実部"と"虚部"である。"実部"の方が左右の耳の間での音の大きさの違いを示すのに対して、"虚部"の方は左右の耳の間での位相差を示すものだ。つまり、ある周波数の音が左右の耳の間でどのような位相差を示すものか、測定しようとするものである。

　左右の耳に対する音の位相差というものは、立体音感を考える上では避けては通れないのだろう。しかし、位相差を処理しようとすると、どうしたらいいものかかなり迷う部分がある。また、今回のようなFFT処理をかけたときに得られる位相を用いて良いものかどうかもよくわからない。というわけで、今回は位相解析処理は後回し、ということにした。

■「ボケ画像」の計算方法について　

　「奥行値(Zバッファ)に従ってボケ量を変化させ被写界深度ブラーを作る場合、0から255の各値に相当するボケ画像(PSF)をFFT→IFFTで作成・合成しようとすると、255枚ものPSF画像とそれぞれを合成するためのマスク画像など、計算コストが掛かりそうですが…」というぐるぐる検索。たいしたことはわかりませんが、ビールを飲みながら適当に書いてみます。あまりに当たり前のことしか書けませんが…。
　計算方法は…、ボケ画像のような畳み込み演算でFFTを使うと(N*Nピクセルの画像であれば)計算量が2N*Nから2N*ln(N)に減るわけですから、力任せに汎用に計算しようとするならば、(よほどPSF画像が処理対象画像に比べて小さいのでなければ)何の迷いもなくFFTを使うのが良いんでしょう(多分)。ただし、何らかの性質がPSF画像やZバッファにあるならば、その性質を利用するのが「手」かもしれません。例えば、非常にPSF画像が非常に(本当に非常に)粗であるとか、Zバッファの値は比較的近くの画素間で連続だとか…。
　また、「255枚ものPSF画像とそれぞれを合成するためのマスク画像」というのはおそらくそこまで必要ではないと思います。もしも、PSF画像を(1+2n)*(1+2n)の正方領域でnが1〜128まで(つまり256*256)まででも、たかだか全部で1/3(1+n)(1+2n)(3+2n)ピクセル程度のサイズですし…。全部前もって作っておいて、それを保持していても大した容量は必要ないですから…。