はこにわオーディオ工学研究分科会 (旧名: バスレフ研究所)

DFT解析シミュレーション

皆様　お正月はどのようにお過ごしでしょうか？

私は年末からMCAP-CR型のシミュレーションプログラミングを楽しんでおり、何とかDFT（離散フーリエ変換）解析による音圧算定の目処が立ってきました。
昨年末に苦しんだC言語での配列のポインタを使ったプログラム変更もウェブの親切な情報のお陰で何とかなりました。また、DFT計算は遅いものの、記述を簡略化出来たので、プログラムが簡潔になりました。高速化しないDFT計算はFFTに比べて遅いとは云え、3200ポイント程度なら速くないノートPCでも一呼吸する間に終わるので実用上は問題はないと思います。6400ポイントでも十数秒で計算が終わります。低域だけなので、これでもオーバースペックなくらいです。
結果を公開する前にもう少し検証が必要ですが、プログラムの内容が非常に簡単なので、DFT部分だけはソースでの公開も検討しています。

音の良し悪しの議論から離れて計算を楽しむというのも良いですね。

謹賀新年

新年明けましておめでとう御座います。

多自由度バスレフ型の日記はブログに移行しました。
今年の重点目標は、多自由度バスレフ型の解析を完成させることです。

今後もよろしくお願い致します。

BSで放送されたワルキューレ

一昨日にミラノ・スカラ座での『ワルキューレ』の公演がNHK BSで放送された。
その時間にテレビにかじり付いて見るわけにはいかないので、録画して視聴した。

オーケストラピットは自分が今までに見たヨーロッパの劇場より奥行きが深い感じがした。そこに、団員がびっしり詰まっての演奏で、張り詰めた緊張感が伝わった。

この劇場でワーグナーというのはあまりピンと来ない感じがしていたが、素晴らしい舞台、素晴らしい熱演、やっぱりミラノ・スカラ座とバレンボイムは凄かった。

決め付けはいけませんね。

プログラミングの目的

昨日はDFTで計算することに成功したのだが、計算量を増やそうとしたら別な問題に当たってしまった。
自分のDFTプログラムでは、初歩のFORTRANプログラムのように、入力データの記憶用に配列を使っていた。
計算結果確認のような小規模のプログラムであれば、配列を使用して問題なかったのだが、配列を増やしていくとエラーになる。これは、想定の範囲内ではあったが、まさかこんなに小さな配列しか使えないとは思っていなかった。1次元配列なので、数万の要素程度はどうってことないと思っていたのだが、1000でも無理だった。

メモリの使用を減らし、大きなデータを扱うには、都度ファイルから読んでいく作戦もあるのだが、実行速度が極端に下がってしまう。そこで、配列を使い、RAMに記憶させて、実行速度を稼ぐという安易な作戦だった。

こういう問題に当たってウェブを検索すると、2色のグループに分類される。『配列の上限値はいくつですか？』というような疑問の意味を理解する人と理解しない人に分かれる。

最大派は、プログラミングそのものが自分のメシの種であり、目的のような人。こういう人は、上記の疑問の意味が分からない。『何故そんなに大きな配列が必要なの？』という疑問を逆に持つ。

自分を含めた少数派は、プログラミングの目的が、計算という人である。物理の問題を数値的に解くには対象を分割する必要がある。有限要素法、境界要素法、差分法などいろいろな手法があるが、こういう問題を解くには、それぞれの要素を配列に格納してしまえば簡単である。だから、配列を安易に扱えれば大変重宝する。コンピュータ内部でどうやってメモリにアクセスするとか、そんなことはどうでもいいのである。例えば流体の挙動を解析するなら、必要なメモリは買ってくる等の投資は惜しまない。自分が決めた手順で自動的に計算してくれ．．．という感じだろう。

最近は、コンピュータの仕組みを熟知した後者グループも居るようなので、自分のようなのは、時代遅れのオールドスタイルになっている。無知で恥ずかしいのだが、もっと安易なプログラムは出来ないのかな？

離散フーリエ変換

標準型MCAP-CRの解析をすすめている。
正弦波入力に対する応答を求め、その後、ランダム波入力に対する応答の計算まで出来たところで、それを周波数分析（FFT分析）しようとしたところで詰まっていた。

GPLライセンスのFFTWサブルーチンを利用したいと考え試行錯誤していた。

最初にFFTWを使える状態にするのに悩んだら、結局OpenSUSEのソフトウェア管理機能を使ってインストールできることが分かった。何だ．．．ここまで1週間

次に、入力データの読込ませ方法が分からない。C言語の基本が分かっていないので、メモリ管理方法から勉強し直して結局何だか分かったような分からないような．．．
直接的に解説してあるページを発見したが、Visual C++を使っているのでそのままは使えない。悩んだが結局何とか計算実行出きるようになった。更に1週間。

計算結果を眺めていても解釈に詰まる。計算結果が正しいか正しくないか分からない。
悩んだ末に、定義通りやってみようかと思い立った。

FFTは面倒だが、DFT（離散フーリエ変換）の定義式は大したことないことに気付いた。自分の研究程度の分量なら実行速度の遅いDFTでも十分か？と思い、一からプログラムを書いてみた。たったの2時間で完成。FFTWの実行結果と十分に一致していた（細かい桁には多少の違いはあるが）。計算は一瞬で終わった。この場合に限れば、実質的にはFFTWを使うより速い。中身が分かっているので、計算結果の出力方法も思い通りだ。

しかし、各ポイントには、正の値だけでなく負の値もある。また、当然ながら、定義通り虚数項もある。
さて、問題は、どうやってこの結果を周波数特性にするのか？？？