今日のおじさん、なに食べました? (仮)

妻の料理と、おじさんの毎日の記録です。ほんのり工学テイスト。

omocha201612.jpg
↑↑↑画像クリックで、我が家のおもちゃを紹介します!↑↑↑
エクセルで周波数解析,FFT(高速フーリエ変換)版
【今日の料理】 2011/1/23 夕食
 今日は,妻と娘(11.4ヶ月)は,午前中は子育て支援センターへ,午後は育児仲間のところへ遊びに行きました.その間,私は,自由時間です.外でランチを食べ,スポーツジムに行き,独身気分を満喫しました.
CIMG7884.jpg
★米飯

★汁物
 昨晩の残り.具沢山の汁物です.

★豚肉味噌漬け
 「サイボクハム」で買った,「豚もも肉の味噌漬け」です.3切れで1200円の,高級品です.大ぶりのもも肉が,味噌で覆われた状態で,販売されています.ホームページ[1]の説明の通り,味噌を大まかに落とし,両面を焼きました.
 とても軟らかく,噛むほどに味が出てきます.ハムに似た食感でした.ただ,私には,ちょっと塩辛すぎました.この味付けならば,半切れでも十分かもしれません.でも,贅沢な味わいで,また買いたくなる品です.
 残った味噌は,もう一度肉を漬けることができるそうなので,試してみたいと思います.
[1]サイボクハム;豚肉みそ漬
http://www.saiboku.jp/item/517.html

★野菜炒め
 ニンニク・もやし・ほうれん草・人参の炒め物です.油は,ごく少量で,蒸し野菜風です.仕上げにバターを少し加えました.
 食べるまでに時間が経ってしまい,冷えてしまったので,いまひとつでした.バターは不要だったかもしれません.

「今日の料理」へのコメント,お待ちしております.
→こちら


↓↓↓料理&育児の応援,お願いします!ポチッと↓↓↓
にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村


【今日の料理工学】 エクセルで周波数解析,FFT(高速フーリエ変換)版
 少し前に,「エクセルで動くフーリエ変換ツール」を紹介しました(→こちら).このツールは,wavファイルを読み込んで,波形を表示して,フーリエ変換を行う,というものでした.いろいろな音の,周波数解析を行うことができます.

 このツールでは,離散フーリエ変換(DFT;Discrete Fourier Transform)を用いています.このため,計算に時間がかかるという弱点がありました.この点が,当時の心残りでした.
 今回,高速フーリエ変換(FFT;Fast Fourier Transform)を用いたツールが完成したので,新たに紹介します.


★エクセルで動くFFTツール
 ツールの画面は,図1のようなものです.エクセル上で動くツールで,マクロ(VBA)を使っています.
<図1>
z20120123z1.jpg

 このツールでは,次のようなことができます.
 ・WAVファイルの波形のグラフ化
 ・波形のフーリエ変換(周波数解析)
 以前のツールとの違いは,フーリエ変換の計算アルゴリズムに,高速フーリエ変換(FFT)を取り入れた点です.データ点数にもよりますが,以前に比べて,10倍程度,計算速度が向上します.

 このツールの利用例は,例えば,こちらの「家電の操作音」編を,ご覧下さい.例えば、次のようなことができます。興味を感じましたら、是非お試しください。
 ・キッチンタイマーや、家電の電子音の周波数を調べる→こちら
 ・ハーモニカのトレモロの秘密(うなり)を調べる→こちら
 ・いろいろな和音を言い当てる、「絶対音感」と勝負!→こちら
 ・自分の「音痴度」を調べてみる→こちら

 また,このツールでは,マクロのソースコードを参照・編集できます.このため,フーリエ変換の計算内容を知ることができますし,自分の好きなように改造もできます.

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
 FFT(高速フーリエ変換)ツールは,こちらです.(xlsファイル)
https://docs.google.com/open?id=0B7vDUGd6_bj-YmE4NzRkMmItY2ViMi00ZTNiLTk0NzMtZDE1YThjOTEyM2Q3
(「Googleドキュメント」が開きます.File→Download Originalでダウンロードできると思います.)
 当ツールは,使用・改変・転載等,全てフリーです.勝手にやって下さい.
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
 当ツールへの改善要望,質問事項などありましたら,お気軽にコメントください.


★使い方 
 使用のためには,マイクロソフトの表計算ソフト「エクセル(Excel)」が必要です.エクセルのマクロのセキュリティを,「中」または「低」に設定する必要があります.(エクセルのメニューバーで,ツール→マクロ→セキュリティ.)

 使い方は,簡単です.数字の書かれたボタンを,順に押していくだけです.
①ファイル指定
 このボタンを押すと,「ファイルを開く」ダイアログが開きます.開きたいWAVファイルを選びます.
 開いたWAVファイルは,左上方にある「再生→」の横のハイパーリンクから開けます.

②WAVデータ読み込み
 続いて,このボタンを押すと,①で選んだファイルのデータが,エクセルに読み込まれます.
 読み込まれたデータは,画面右上に,波形のグラフとして表示されます.

③解析区間切り出し
 ②で表示された波形を見ながら,どの部分をフーリエ解析(周波数解析)したいか決めます.
 そうしたら,このボタン③を押します.
 すると,「切り出し開始時間は?」と聞かれるので,解析開始したい時間を入力します.
 続いて,「切り出し間引き間隔は?」と聞かれます.これは,1のままでOKです.
 そうすると,画面右中央に,解析対象となる波形が現れます.

④FFT/DFT
 最後に,このボタンを押すと,③の範囲がフーリエ変換されます.
 結果が表示されるまで,少し時間がかかります(おおむね数秒~数十秒).
 周波数解析の結果は,画面右下のグラフに表示されます.
 なお,切り出したデータ数が2のべき乗(N=256,512,1024,2048,4096,…など)の場合に限り,FFTとなります.これ以外のときは,DFTで計算されるので,計算が極端に遅くなります.(切り出すデータ数の変更方法は,後述.)


★注意点など
②で,
 あまり長いデータ(50万点超)は,途中までしか読めません.
 また,データが5000点以上の場合は,5000点程度になるよう,間引きながら読み込みます.
 
③で,
 切り出された箇所が,思った箇所でなかった場合,もう一度③を押せば,やり直しできます.
 切り出すデータ点数を増やしたい場合は,「切出点数」のセルの数字(初期値:1024)を変更します.データ数が2のべき乗(N=256,512,1024,2048,4096,…など)の場合に限り,FFTとなります.これ以外のときは,DFTで計算されるので,計算が極端に遅くなります.
 周波数分解能を上げたいときは,「切り出し間引き間隔は?」に,2以上の数字を入れます.そうすると,解析範囲(時間)が長くなって,周波数分解能が上がります.周波数分解能[Hz]=1÷解析範囲[sec],です.
 また,モノラルのみ対応です.(ステレオの場合,どうなるか分かりません.)

④で,
 フーリエ変換には,通常は,適当な「窓関数」を用いるべきのようです.しかし,本ツールでは窓関数は「なし」(矩形窓)を標準とします.「FFT/DFT Hamming」のボタンを押すと,「ハミング窓」での解析ができます.他の窓関数が必要な場合は,マクロを改造してください.


 興味のある方のために,高速フーリエ変換がどのようなものか,以下で説明しておきます.

★離散フーリエ変換(DFT)と高速フーリエ変換(FFT)
 いろいろなフーリエ変換については,過去の記事(→こちら)で,簡単に説明してあります.
 離散フーリエ変換と,高速フーリエ変換は,計算の内容は同じです.しかし,計算アルゴリズムが違います.FFTでは,解析されるデータ数を2のべき乗に制限する代わりに,DFTに比べて非常に短時間で計算を完了できます(例えばデータ数N=1024個のとき,必要な掛け算の回数はDFT:約100万回,FFT:約1万回).

 以下では,離散フーリエ変換から始めて,高速フーリエ変換の計算手法を導きます.

 まず,離散関数g(k)をsinとcosで表す「離散フーリエ逆変換」は,次式1です.Nは,データ数です.また,k=0,1,…,N-1,です.
<式1>
20120123s01.jpg
※多くの解説書では,e(自然対数の底)のj乗(複素数)の形で,フーリエ変換を表しているようです.しかしここでは,「とっつきやすさ」を重視して,sinとcosを用います.そのかわり,計算は少し煩雑になります.

 式1の係数A0,An,Bnを計算する式が,次式2の「離散フーリエ変換」です.
<式2>
20120123s02.jpg

 高速フーリエ変換では,データ数Nが,2のべき乗のときに限定します.すなわち,N=2^rです(r:整数).そうすると,データ数Nは,必ず偶数個です.
 式2を,次式3のように変形します.ここで,後の計算での混乱を避けるため,式2のAnをAnGと表しています.
<式3>
20120123s03.jpg
 式3では,離散関数g(k)の,偶数番目をe(k),奇数番目をh(k)と,新たに定義しています.e(k)およびh(k)のデータ数は,それぞれN/2個になります.

 ここで,e(k)とh(k)の離散フーリエ変換を考えます.式1でg(k)=e(k)またはh(k),N=N/2とおくと,次式4が得られます.ここで,e(k)およびh(k)のA0,An,Bnの組を,A0E,AnE,BnEおよびA0H,AnH,BnH,と表します.
<式4>
20120123s04.jpg

 式3と式4から,n=0,1≦n≦N/2-1については,次式5が導かれます.
<式5>
20120123s05.jpg
 n≧N/2の場合は,式5を使おうとすると,AnE,BnE,AnH,BnHを計算できません.これらの係数は,1≦n≦N/2-1の範囲でしか,定義されないからです.そこで,別途,算出方法を考える必要があります.
 まず,n=N/2については,式3でn=N/2とおいて,次式6を得ます.
<式6>
20120123s06.jpg
 次に,N/2+1≦n≦N-1については,式3でj=n-N/2とおくことで,次式7を得ます.
<式7>
20120123s07a.jpg
20120123s07b.jpg

 以上の式5,6,7から,An0,AnG,BnGは,次式8で計算できます.
<式8>
20120123s08.jpg
 すなわち,g(k)の離散フーリエ変換は,e(k)とh(k)の離散フーリエ変換で,表すことができます.これが,高速フーリエ変換の,基本的な考え方です.
 

★高速フーリエ変換(FFT)が高速な理由は
 離散フーリエ変換の式は式2ですので,係数An,Bnを定めるために,概ね2×N^2回のかけ算が必要です(g(k)とcos,sinのかけ算のみに注目しています).したがって,g(k)およびe(k),h(k)の係数An,Bnの算出に必要なかけ算の回数は,次のようになります.
 ・g(k):2×N^2 回
 ・e(k):2×(N/2)^2=2×(N^2)/4 回
 ・h(k):2×(N/2)^2=2×(N^2)/4 回
・e(k)とh(k)の両方を計算:2×2×(N^2)/4=2×(N^2)/2 回
 すなわち,式8のように,g(k)をe(k)とh(k)に分解することで,かけ算の回数を,2×(N^2)回から,2×(N^2)/2回に,半減できます.(実際には,式8の中にもかけ算がありますので,ここまでは減りませんが,大幅に減ることは変わりません.)

 コンピュータの計算の速度は,概ね「かけ算」の回数によって決まるようです.したがって,かけ算を減らせれば,計算スピードを向上できます. 

 ただ,式8を使うだけでは,かけ算の回数を半分にできるに過ぎません.そこで,式8のe(k)を新たなg(k)として,これを新たなe(k)とh(k)に分解します.そうすると,式8のe(k)も,半分のかけ算回数で算出できます.このように,関数をどんどん分割していくことで,かけ算回数をどんどん減らしていけます.最後には,データ数が1個のg(k)が残ります.このような分割を続けていくために,データ数Nが,2のべき乗である必要があるのです.

 式2から,N=0のときは,g(k)の離散フーリエ変換は,次式9となります.
<式9>
20120123s09.jpg
 結局,高速フーリエ変換には,式8と式9があればよいことになります.これが,高速フーリエ変換の計算手順です.離散フーリエ変換の式2は,直接は使用しません.しかし,計算している内容は,式2と同じです.この計算式のプログラム化については,上のツールのソースコードを参照ください.(式8の計算ルーチンを,再帰的に呼び出すことで,実現できます.)

 高速フーリエ変換(FFT)と離散フーリエ変換(DFT)のかけ算の回数の比は,より正確には,次のようになるそうです.
 FFT:DFT=(Nlog2N):(N^2)=log2N:N
 具体的なデータ数Nを用いると,次のようになります.
 ・N= 256(=2^8 )のとき:FFTのかけ算回数≒DFTの1/30
・N=1024(=2^10)のとき:FFTのかけ算回数≒DFTの1/100
 ・N=4096(=2^12)のとき:FFTのかけ算回数≒DFTの1/300
 データ数が多いほど,FFTの高速化の恩恵が,大きくなります.

 ちなみに,今回紹介したFFTツールで試してみると,N=4096のとき,計算の所用時間(計算ボタンを押してから,結果が表示されるまでの時間)は,次の通りでした.
 ・DFTで計算:約38秒
 ・FFTで計算:約 3秒
 FFTだと,約10倍の速度になっています.上述の理屈では,約300倍のはずですが,実際には,データ読み込みと結果の表示に結構な時間を使っているようで,ここまでの高速化は実現できないようです.


 当ツールへの改善要望,質問事項などありましたら,お気軽にコメントください.


【今回の結論】
 エクセルで動く,FFT(高速フーリエ変換)ツールを紹介しました.
 FFTは,DFT(離散フーリエ変換)と同じ計算を行いますが,計算手順が違います.大幅な高速化が可能です.


今回は,完全に趣味の話です.申し訳ありません.でも,これだけやらないと,心残りだったのです.
応援くださると嬉しいです.
↓↓↓↓↓↓↓↓
にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

関連記事
スポンサーサイト
omocha201612.jpg
↑↑↑画像クリックで、我が家のおもちゃを紹介します!↑↑↑









管理者にだけ表示を許可する


こんばんは。
趣味の観点から、さっそく周波数解析ツールを試してみました。文系なので(この言葉は禁句だったか?)数式の吟味などブラックボックスの中は完全スルーで、いきなり音楽ファイル(しかもステレオ)を読みこませました。カセットテープ保存のアナログデータをデジタル化したものです。
データ数が多い割には2分ほどで計算してグラフを出してくれました。しかし残念なことに周波数は0~2250Hzの低く狭い範囲にしか分布していないという結果が出ました。これはさすがに正解とは思えません。おそらくデータの入力時に不備があったのでしょう、ステレオデータであるとか、収録時間であるとか。
今度はアナログにして試してみます。
趣味に使えそうなプログラム、ありがとうございました。
niwatadumi | URL | 2012/01/24/Tue 23:03 [編集]
> niwatadumi さん
 こんばんは.
 早速使ってくださり,ありがとうございます.でも,うまく動作しなかったようで,残念です.
 計算に2分かかったということは,相当データ数が多いような気がします.私が確認したのは,せいぜい録音時間が10秒くらいのデータでした.今回の場合,プログラムで想定したデータ数の上限を超えてしまったのかもしれません.(データが2250Hzで切れているとなると,この可能性が高そうです.)
 懲りずに使ってみて下さると,ありがたいです.
マツジョン | URL | 2012/01/26/Thu 00:18 [編集]
こんばんは。
お邪魔致します。
私の様な者が使いこなせるかは分かりませんが、とても興味深いと思いコメントさせて頂きました。
ニコニコ | URL | 2012/01/29/Sun 23:42 [編集]
> ニコニコ さん
 こんばんは.
 訪問ありがとうございます.つたないツールですが,使ってみて下さると,とても嬉しいです.
マツジョン | URL | 2012/01/30/Mon 23:20 [編集]
失礼します。
質問なのですが、上図1においてFFTを行う前の二つのグラフで、
横軸の時間[ms]に対して縦軸は何を表しているのでしょうか。(-2から2までの目盛)
またその縦軸に対応してプロットされているデータとは別に、nsamplの隣にプロットされているデータが何を示しているのかを教えていただけるとありがたいです。
素人質問ですみません
| URL | 2016/02/14/Sun 00:29 [編集]
こんばんは。ご質問ありがとうございます。

 縦軸は、信号の強さです。

 まず、nsampleの横のデータ(3列目)が、wavファイルの信号の生データ(ビット表示)です。8ビットの音声データの場合、0から255の数字で示されます。この数字は符号付き8ビットなので、信号としては-127から+128を表しています。

 そして、グラフのプロットに使っているデータ(7列目)が、この生データを-1~+1の範囲に加工したものです。具体的には、(3列目のデータ-128)/128で計算しています。±1が、録音レベルの最大値となります。

以上、お役に立てましたら幸いです。
マツジョン | URL | 2016/02/16/Tue 21:47 [編集]
信号レベルが-1近辺になる
ダウンロードして使ってみました。
それらしく動作している様ですが、何故かどのWavを読み込んでもすべて-1を中心として上下に変化する波形として取り込まれています。
マクロのOpenFilesの中の

Sub WriteData_all(wdata0() As Integer, fs As Long, n_data As Long, maxlevel As Long)



Cells(R_, C_ + 3) = wdata0(i) / maxlevel * 2 - 1 '*** データ(最大±1)

を、

Cells(R_, C_ + 3) = wdata0(i) / maxlevel * 2 '- 1 '*** データ(最大±1)

とし、

Sub WriteData_trim(wdata0() As Integer, fs As Long, n_data As Long, maxlevel As Long)

も、同様に -1 をコメントアウトすればデータ波形は0を中心に取り込まれます。

ここで-1としているのはどういう意味が有るのでしょうか?
また、この-1をやめた場合、どういう問題が有るのでしょうか?

Wavの構造及びFFTの処理を完全に理解出来ていないので、おかしな質問かもしれませんが、よろしくお願いします。
booken | URL | 2017/05/10/Wed 23:31 [編集]
> booken さん
 ご質問ありがとうございます。

 wavの元データは8ビットの数字で、0~255の範囲で表されます。これを128で割って1を引いて、-1~+1の範囲のデータに変換する(スケールを変える)のが、ご指摘の数式の意図です。
 ただ、PC・エクセル・wavファイルのバージョン・環境・設定等によっては、0~255を、最初から符号付きの-127~+128として読み込まれることも多いようです。この場合、128で割るだけで、-1~+1の範囲のデータに変換され、1を引く必要がありません。bookenさん指摘の通り、-1をコメントアウトするのが妥当です。

 なお、-1してもしなくても、FFT結果には実質的に影響はありません(直流成分、すなわちf=0Hzの成分が変化するだけ)。wavグラフ表示の不自然さだけの問題となりますので、-1するかしないかは、お好みでお選びくだされば、と思います。

ご使用くださいまして、ありがとうございます。
マツジョン | URL | 2017/05/12/Fri 06:41 [編集]
ご回答ありがとうございました。
WAVのデータ構造が符号なし、VBAではInteger(符号付き)で読み込んでいるが、符号なしで読み込む場合と符号付きで読み込む場合があり、符号なしで読み込んだ場合は-1が必要だが、符号付きで読み込んだ場合は-1不要という事と理解しました。
私の環境では符号付きで読み込むようなので、-1をコメントアウトして利用させていただきます。
ありがとうございました。
booken | URL | 2017/05/12/Fri 22:40 [編集]
> booken さん
 解決したようで嬉しいです。ご活用くださいませ。
マツジョン | URL | 2017/05/20/Sat 10:29 [編集]



トラックバック
TB*URL







Copyright © 今日のおじさん、なに食べました? (仮). all rights reserved.