あちゃぴーの自転車通勤

空気中の音速は「秒速 340m」と覚えておいても大抵問題ないのだが、空気の温度によって音速は大きく変化する。また風向きとかでも変化してしまう。この 340m/s(秒)は摂氏 15 度のときの速さになる。温度で速度が変わるのは、空気を構成する分子の運動速度が温度が上がるほど速くなるから。 空気は温度以外に、湿度や気圧などの変化もある。音速は気圧の変化には影響を受けないが、音圧は影響を受ける。また湿度にも影響を受けるが、温度変化に比べれば、影響は少ないため通常は考慮しない。 温度が変化すると密度も変化する。密度が変化することで音速が変わるという言い方もできそうだが、なかなか手ごわい話になりそうなので、ここでは触れない。下のプログラムでは、温度から密度を計算している。 音速の計算 音速は波動方程式から導いているが、細かくやると大変なので、流体中の音速の方程式からスタートして、空気専用にまで式を変形させてみたい。 まず元となる式はこれで、結構シンプル。 c 音速は通常これを使う 単位は m/s K 体積弾性率 ρ 密度 体積弾性率を空気の場合に限定してみる。 γ は比熱比 P 0 はヘクトパスカル 標準大気圧 1013.25 hPa となるが、単位をそろえるために 101325 Pa(N/m 2 ) これを使う。 γ 比熱比を分解してみる γ = c p /c v c p 低圧比熱 1.006 c v 定積比熱 0.717 1.006 / 0.717 = 1.403068340306834 今度は温度によって、密度 ρが変化するので、それに対応させる。温度と気圧から求める。 273.15 は摂氏温度(セルシウス温度)を絶対温度変換するため。 t は摂氏温度 P 0 は気圧 101325 Pa(N/m2) R は気体定数 2.87 これで温度ごとに音速が計算できる。上記をすべてひとつの式に入れるとこうなる。整理していないので、かなり汚い。中身で何を計算しているのか分かりやすいので、下の計算プログラムは、この式をそのまま使っている。 実は式を整理していくと下のように気圧は消えてしまう。音速は気圧に影響を受けないことを物語っている。 また温度以

「中学生ぐらいから学ぶフーリエ変換」ということで、中学生ぐらいでも分かるようにまとめようとは思っているが、まずは、書きたいことを書いて、後々改善していきます。分かりやすく書くことほど難しいものはない。 今回はいきなり、縦波(たてなみ)と横波(よこなみ)の関係について。振動は空気や物体などの媒質(medium)を伝わっていくのだが、振動は縦波(longitudinal waves)と横波(transverse waves)に分類できる。人が音として認識しているのは縦波となる。横波は、ねじっても戻る(せん断応力のある)ような固体などは伝わることができるが、気体や液体ではほとんど無視できるもの。一般に音として扱うのは縦波となる。また縦波は疎密波とも呼ばれる。 空気中を伝わる音(縦波)とは、実際には空気の伸び縮みで、細かく見ていくと空気の粒子が振動方向に対して細かく往復運動して、次々に伝えて波となる現象。ひとつの粒子が遠くまで移動しているわけではない。この縦波を人間が工学的に音を扱う場合、縦波のままだと不便なため、通常は横波に変換して扱っている。 下図(canvas &amp JavaScript)はサイン波を縦波と横波で表現している。スライドバーで wave length (周期)と gain (音量)をいじることができる。(HTML5に対応していない古いブラウザでは表示されないか、妙な表示がされているかもしれない。) グリーン：サイン波を離散(ある等間隔時間の音量だけを見ている)信号にしたもの。 レッド：信号を元に縦波に変換したもの。密度の高いところ低いところがあるのがわかる。これはそのまま音圧(sound pressure)となっている。この音圧の繰り返しが音として認識される。下の表示に限っては、グリーンの棒の長さを 1/10 に縮小して、それを時間軸に置き換えている。信号が 0 のところは移動量 0 であり、信号がプラスであれば、プラス時間側へ移動。マイナスであれば、マイナス時間側に移動している。 ブルー：音圧をサイン波で表現したもの。密度の高いところがプラスの最大値で、密度の低いところがマイナスの最低値。信号のサイン波とは 1/2π ずれている。 縦波と横波の関係 wave length

先日のメトロノーム にボリュームを加えてみた。コードも少しシンプルに変更。自分で使う分には、これで十分だわ。exe のサイズが 9KB というところと、省メモリ、起動の速さがポイント。それにしても GUI がないアプリはソースコードが短くなって楽。GUI のダラダラしたソースコードは公開するのも気が引けるけど、これぐらいなら見通しもよいと思う。 コンパイルされた exe ファイルをダブルクリックで起動。初期値 120 のテンポで 0.9 のボリュームで自動再生。テンポを変更するときは数値を入力して Enter する。ボリュームを変えるときは v0.5 のように打って Enter する。終了するときは q Enter とする。 OpenAL metronome2 volume ソースコード /* metronome2.c * February 17th, 2013 * Namagi Products メトロノーム ボリューム付 * compile win * gcc metronome2.c -lOpenAL32 -lalut -s -o metronome2 * compile linux * gcc metronome2.c -lalut -s -o metronome2 * BPM（Beats Per Minute）60bpm */ #include &ltstdio.h&gt #include &ltstdlib.h&gt #include &ltmath.h&gt #include &ltstring.h&gt #include &ltAL/alut.h&gt void input(double *amp, double *bpm){ char str[32]; int i; /* 入力 */ for(i=0; i&lt3; i++){ printf("input:\t"); fgets(str, 32, stdin); if(str[0]=='q'){ *bpm = -1; break; /* ボリュームの調整 */

OpenALで音が出せたので、メトロノームを試しに作ってみる。出来たアプリケーションは 9KB と小さく起動は高速。ダブルクリックすると、コマンドプロンプトが開いて初期値 120 のテンポで鳴り続ける。他のテンポにしたいときは、テンポを入力して Enter する。小数でも入れられるようにした。テンポの範囲は 20 から 300 まで。終了したいときは q を入力して Enter。メトロノームの音はプログラムの中でサイン波を作ってみた。音にパワーがあまり感じられないので、今後メトロノームに合った波形を作ってみたいと思う。 翌日ボリューム調整可能なバージョンも作ってみた。 OpenAL metronome ソースコード /* metronome.c * February 16th, 2013 * Namagi Products メトロノーム * gcc metronome.c -lOpenAL32 -lalut -s -o metronome * BPM（Beats Per Minute） */ #include &ltstdio.h&gt #include &ltstdlib.h&gt #include &ltmath.h&gt #include &ltAL/alut.h&gt double input(void){ char str[32]; int i; double bpm = 0; for(i=0; i&lt3; i++){ /* tempoの入力 */ printf("input\tBPM:\t"); /* 整数入力 */ fgets(str, 32, stdin); /* 数値判定 数字でなければ再入力 3回まで */ if(str[0]=='q'){ return -1; }else if(isdigit(str[0]) != 0){ bpm = atof(str); break; } } /* Tempoの範囲調整 */ if(bpm &lt 2

10週目にもなると、あえて書くこともなくなりつつある。ギターを始めて一気にうまくなるのは2ヶ月目ぐらいまで。その後は地味な上達という感じ。今現在で、コードを見れば大抵のコードは押さえられるし、コードネームからコードを見つけて押さえるということはできるようになった。新しい曲をやるときでも、コードさえあれば、数分後にはなんとなく弾けている。細かな部分では改善の余地は当然あるのだが、練習の内容を文字にして残しておくほどでもないと思えた。 作曲 10週目以降は、ほとんど練習らしいことはせず、なんとなく作曲をやっている。コード進行を作って、その上にメロディを乗せていくという手順。無数にある組合せの決定が難しい。イメージがどれだけしっかりあるかがポイントで、それが判断基準となる。1曲目がまとまるのはいつのことやら。 歌 弾き語りをする上では歌は必要。これが結構難しい。歌がうまければ、ギター下手でも成立するけど、逆は残念な感じ。ギターと同じように練習というか、トレーニングが必要かもしれない。 9週目＜     ＞7ヶ月振りにギターを再開してみる ギター練習

ANSI C89 の範囲から外れてしまうが、音をリアルタイムに扱ったり、ゲームのちょっとした時間調整用に Sleep関数 を使ってみる。 UNIX、Linux では標準的に sleep関数 を使うのだが、Windows 上では使えないことが判明。代替として &ltwindows.h&gt にある Sleep関数 を使う必要がある。注意点は単位の違い。標準は sleep(sec) だが、Windows は Sleep(msec) となっている。 windows.h Sleep を使って一時停止 #include &ltstdio.h&gt /* Win32 API Sleep(msec)を使うために必要 */ #include &ltwindows.h&gt /* 標準sleep(sec)をWinのSleep(msec)に置き換え */ #define sleep(n) Sleep(n * 1000) int main(void){ int i; for(i=5; i&gt=0; i--){ printf("%d\n",i); sleep(1); } return 0; } 実行すると、5から1秒ごとにカウントダウンして、0になったら終了する。 5 4 3 2 1 0 C言語 ANSI C89 Meadow & MinGW GCC 目次

OpenALでサイン波を生成＆再生させてみた。結構簡単にできるので面白いかもしれない。 OpenAL サイン波生成＆再生 /* * sine.c * シンプルなサイン波を鳴らす * February 14th, 2013 * Namagi Products * gcc sine.c -lOpenAL32 -lalut -s -o sine */ #include &ltstdio.h&gt #include &ltmath.h&gt #include &ltAL/alut.h&gt int main(int argc, char *argv[]) { alutInit(&argc,argv); ALuint buffer, source; alGenBuffers(1, &buffer); alGenSources(1, &source); int i, hz; /* サンプリング周波数 */ int samplingfreq = 44100; /* 点間のステップ */ float step = 2 * M_PI / samplingfreq; /* 角速度 */ float w = 0; /* 出力最大 1 */ float amp = 0.5; char str[8]; for (i = 0; i &lt 3; i++) { /* 周波数の入力 */ printf("sine wave frequency? 20-10000Hz\n"); /* 整数入力 */ fgets(str, 8, stdin); /* 数値判定 数字であれば抜ける 数字でなければ再入力 3回まで行う 1文字目しか見ていないけど */ if((isdigit(str[0]) != 0)) { hz = atoi(str); break; } } /* 周波数の範囲調整 */ if (hz &lt 20 ) { hz = 20; } e

C言語 ANSI C89の範囲だけでは直接音を出すアプリケーションは作れそうもない。Windows のライブラリを使えば簡単に実現できそうだったけど、OpenAL(Open Audio Library) というマルチプラットフォームなオーディオ API があるので、それをインストールしてみた。 インストール OpenAL公式ページから以下をダウンロードする。 http://connect.creativelabs.com/openal/default.aspx OpenAL11CoreSDK.exe freealut-1.1.0-bin oalinst.exe OpenAL11CoreSDK.exe ダウンロード先は Creative Labs: Connect > OpenAL > Downloads Windows へのインストールは OpenAL11CoreSDK.exe をダブルクリックして行う。デフォルトのインストール先は C:\Program Files\OpenAL 1.1 SDK となっている。これは好きな場所で構わないと思う。どうも VC++ での開発を意識しているようなので、MinGWを使う場合は、いろいろファイルを移動する必要がある。 C:\Program Files\OpenAL 1.1 SDK\libs\ 内のファイルすべてを C:\MinGW\lib\ へ移動。 C:\Program Files\OpenAL 1.1 SDK\include\ 内のファイルすべてを C:\MinGW\include\AL\ へ移動。(AL フォルダは自分で作成) そのほかの OpenAL11CoreSDK のファイルは使わないので、捨ててしまっても問題ないようだ。 freealut-1.1.0-bin ダウンロード先は Creative Labs: Connect > OpenAL > Downloads > ALUT 解凍し、その中の下記ファイルを移動する。 \lib\ alut.dll は C:\WINDOWS\system32\ へ入れる。 \lib\ alut.lib は C:\MinGW\lib\ へ入れる。 \include\AL\ al

strncmp()を使った文字列の比較。wavヘッダを読み込むときに頻繁に利用するので、改めて基本的な動作を確認。 文字列の比較 #include &ltstdio.h&gt #include &ltstring.h&gt void cmp(char *str1, char *str2) { int bytecount = 4; /* 文字列が等しいならば 0 (byte数を指定可) */ if (strncmp(str1, str2, bytecount) == 0) { printf("The strings are identical." "\"%s\" \"%s\" %d-bytes", str1, str2, bytecount); } else { printf("The strings are different." " \"%s\" \"%s\" %d-bytes", str1, str2, bytecount); } } int main(void) { char str[][8] = {"image", "imagine"}; cmp(str[0], str[1]); return 0; } 頭から4byte分同じであれば、5byte以降は違っても同じ文字列とみなすプログラム。 The strings are identical."image" "imagine" 4-bytes C言語 ANSI C89 Meadow & MinGW GCC 目次

出来れば使わない方がよさそうな関数の再帰呼び出し。 再帰呼び出しは、ある関数が処理中に、自分自身を呼び出す処理のこと。 通常はfor文で実現できるので、コードの見やすさ、処理スピードを考えると、なるべくforで実現した方がよさそう。中には再帰でしか実現しにくいアルゴリズムもあるようだが。 以下は再帰呼び出しをあえて使ってみた。 再帰呼び出し #include &ltstdio.h&gt /* 再帰呼び出し */ void recursivecall(int n) { if (n &lt 5) { printf("recursivecall() enter\t%d\n", n); recursivecall(n + 1); printf("recursivecall() exit\t%d\n", n); } } int main(void) { recursivecall(0); return 0; } ecursivecall()という関数を再帰呼び出しするプログラム。メイン関数からはじめにrecursivecall()が呼ばれてインスタンス化される。さらにrecursivecall()の中でrecursivecall()が呼ばれる。これを4回繰り返したら、ひとつずつ抜け出してくるという仕掛け。再帰呼び出しというと自分自身を呼び出すと書かれている本が多いけど、どうも新規に自分と同じものを作り出して、入れ子になっていると考えた方が自然なような気がする。 recursivecall() enter 0 recursivecall() enter 1 recursivecall() enter 2 recursivecall() enter 3 recursivecall() enter 4 recursivecall() exit 4 recursivecall() exit 3 recursivecall() exit 2 recursivecall() exit 1 recursivecall() exit 0 再帰呼び出しで黄金比を求めてみる #inclu

fgets C言語でコマンドプロンプトから数値を入力する方法はいくつかあるのだが、ここでは fgets を使って文字列として読み込み、数値に変換する方法を試してみる。 fgets 以外にも scanf を使った方法があるが、予期せぬ入力に対応しきれないため、なるべく扱わない方がよいらしい。 atoi, atof 文字列を数値に変換する方法もいくつかるが、ここでは整数に変換する atoi と、小数に変換する atof を扱う。これらを利用するためには、ヘッダファイルの &ltstdlib.h&gt をインクルードする必要がある。 整数の入力 fgetsとatoi コマンドプロンプトから入力した数値(文字列)をintに変換する例。分かりやすく、入力した文字列も表示するようにしてみた。 #include #include int main(void) { char str[32]; int num; printf("整数を入力してください。\n"); /* 整数入力 */ fgets(str,31,stdin); /* charをintへ変換 */ num = atoi(str); printf("入力された整数 int: %d\n", num); printf("入力された文字列: %s\n", str); return 0; } intに入るギリギリの大きさの数値を入れてみる。問題なく変換されているのが分かる。 整数を入力してください。 2147483647 入力された整数 int: 2147483647 入力された文字列: 2147483647 次にintに入りきらない大きさの数値 2147483648 を入れてみると、intに変換された後は違った数値になっているのが分かる。実際に使用する場合は、なんらかの対策が必要となる。 整数を入力してください。 2147483648 入力された整数 int: -2147483648 入力された文字列: 2147483648 次に予期せぬ文字列が入った場合どうなるかを試してみる。123abc456def789 という文字列を入れてみる

改めてWavファイルのヘッダを読み込むプログラムを作ってみる。ここにあるプログラムは完成したものではなく、予備実験のサンプル。関数の配置を試行錯誤しているところ。以前はファイルから直接一字一字読み込んでいたけど、今回は、まるっとメモリに読み込んでしまうことを考えてみた。この後の処理はそのうち書こうと思う。 プログラムの内容 コンパイルして出来たexeファイルにwavファイルをドラッグ＆ドロップすると、ヘッダを1byteごと16進数で表示する。ただし現状では60byte固定。足りない場合は2倍にする仕組みを今後入れていく。またヘッダの解析を入れるとゴチャゴチャするので、今は入れていない。 プログラムの流れは、機能ごとに関数にわけて、main関数は各関数に指示を出すだけにした。これによって柔軟性のある開発が可能になる。共通で使う値などは構造体に管理する。main関数でインスタンスすることで、他の関数からも使えるようにしてみた。構造体に項目を追加して実用になるプログラムにするつもり。 wavヘッダ読込プログラム /* * header_r.c * February 10th, 2013 * Namagi * ヘッダを読むための予備実験 ヘッダをそのままメモリに格納 */ #include &ltstdio.h&gt #include &ltstdlib.h&gt /* バッファ構造体 */ typedef struct { /* ファイルポインタ */ FILE *fpr; /* ヘッダを丸々入れておく そのポインタ */ unsigned char *buff_header; } BUFF; /* ヘッダ情報構造体 */ typedef struct { /* ヘッドのサイズを入れる */ int headsize; /* ファイルのパスを格納 */ char *file_path; } HEAD; /* ヘッダリーダー 構造体に入れていく処理 */ void header(char **argv, BUFF *buff, HEAD *head) { int i; /* 初期値 足りなかったら倍にする そのうち実装 */ head-&gthe

gcc のバージョンを確認 まずは現在のバージョンを確認する。コマンドプロンプトや Meadow の shell などから gcc -v と打つとバージョンが表示される。 現在の 4.7.0 が確認できる。 gcc の最新バージョンをダウンロード MinGWで使えるgccの最新バージョンは 4.7.2 なので、これを下記URLからダウンロードする。 http://sourceforge.net/projects/mingw/files/ ページを開くと、このような画面が出る。MinGW をクリックする。 さらにgccの最新版の gcc-4.7.2-1 まで辿っていく。 Home / MinGW / Base / gcc / Version4 / gcc-4.7.2-1 この階層にバージョンアップに必要なファイルがリストアップされている。その中の gcc-core-4.7.2-1-mingw32-bin.tar.lzma 2012-11-05 9.9 MB をダウンロードする。 gcc の最新バージョンのインストール ダウンロードしたファイルを解凍すると以下の3つのフォルダができる。 bin, lib, libexec それぞれをMinGWのインストールされているディレクトリに上書きして完了。 再度コマンドプロンプトから gcc -v と打ってバージョンを確認してみる。 最新の gcc version 4.7.2 (GCC) が確認できたので成功。 C言語 ANSI C89 Meadow & MinGW GCC 目次はこちら

y = 2x こんな式を中学で習う。x に何かを入れると y が求められるという。 一般的に入力は x で、出力は y と決められている。逆はだめ。何で？ 求まるジャン！ 高校になると y → f(x) このように置き換えましょう。と言われてしまう。y は何だったの？という疑問が出てしまう。とりあえず言われたとおり、上の式に当てはめるとこうなる。 f(x) = 2x f は function の略で単なる名前に過ぎない。必ずしも f(x) である必要はなく、g(z) でも h(t) でも構わない。 専門書は、この記述で書かれているので、これを使うと、ちょっと大人っぽい気になったりもした。 でも改めて考えてみると不思議な気もする。わざわざこんな書き方しなくてもいいじゃんと。 強い疑問でもなかったので長い間放置状態が続く。 最近C言語をやってみると、f(x) に対するモヤモヤとしたものが取れ始めた。 C言語では関数単位でプログラムを作る。まさに数学の関数と同じ。書き方もそっくり。例えば、 int function(){} function が関数名で ( ) の中には引数が入る。つまり f(x) の x が引数に当たる。数学と全く同じ構造になっている。そして処理内容は { } の中に書く。C言語ではイコールは使われないが { } が同じ働きをする。「f(x) = 2x」をCで書くと下のようになる。 int function(int x) { return 2 * x; } 「function()」という名前の関数に「x」という整数を入れると、「2を掛けて戻す」という内容。 C言語と同じ感覚で数式を見ると、数式の意図が明確に見えてくる。確かに y を使っていたら、y に入れて x を導こうなんて考えてもおかしくないが、f(x) だったら、その気もなくなるというわけ。今度は何で「y」を使っているの？という疑問が出てしまった。パラドックスか？ 予想するにはグラフとのかかわりだろう。 中学生ぐらいから学ぶフーリエ変換 目次はこちら

はじめに 音処理について、数学的に学習して行こうというコーナー。また実際の処理としてプログラミングも積極的に利用するつもり。以前から書こうと思っていた内容なのだが、ようやくスタート。タイトルは、なんとも中途半端だが、その内容は息子(小学生)向けの補足資料的な扱いが強い。内容的には中学生でも理解できるようにしたいが、実際には高校数学の範囲もやることになる。ということであまり範囲を決め付けず、音処理に必要な数学ならば大学レベルでも学習したいとは思う。難しいと思えたことは詳しく解説するつもり。たぶん参考書の解説とは方向が違うので、大学の先生には怒られそうな内容になるかもしれない。調べまくった正しい解説ではなくて、こう思って使っているレベルの話が多くなるはず。ひょっとしたら誤解を招くかもしれない。 数学そのものの学習という意味合いも少しはある。音処理という方向を定めておくことで、数学やプログラミングの目的がはっきりするので、抽象的な数学でも具体的な応用例が実感できて分かりやすいと思われる。目標としてはフーリエ変換ぐらいまでとしたい。これを理解するには高校数学の広範囲の知識が必要だし、さらにプラスアルファも必要だと思う。それなりの量という感じだ。 音そのものについてや、基礎的な数学からスタートするので、何年かかるか分からないし、途中でやめてしまうかもしれない。とりあえず実際に子供に教えつつ更新して行くので、ボチボチとやっていきます。 数学とプログラミング 「コンピューター」は別の言い方をすると「計算機」と言う。パソコンやスマートフォンが日常的に使われている現在は、計算機というイメージは失われているかもしれない。それでもノイマン型コンピューターの中身は2進数の計算しかしていない。 プログラミングを実行するにはコンピューターが必要。コンピューターは計算機なのだから、プログラミングは計算式となる。計算式は算数や数学の数式とそれほど変わらない。実際にプログラミングを組めば、数学とプログラミングがとても近い関係だということが分かると思う。数学が理論的なことを教える座学だとすれば、プログラミングは数学を利用した実習となる。見方を変えれば数学は抽象的な世界であり、プログラミングは具体的な世界とも言える。この違いは膨大な計算をするか

FIRフィルタを作る場合には必須の手法。サンプル数が長い、もしくは不明な入力信号と、サンプル数が分かっている(少ない)フィルタをFFTで畳み込んでみる。手法は overlap-add method を使う。これは入力信号を一定のサンプル数で切ってから、FFTにかけて変換し、複素数としてフィルタと演算後、iFFTで戻す。戻った細切れの信号はフィルタ分-1長くなるので、その部分をオーバーラップしながら足して、出力信号とするもの。 下図は、この記事のプログラムを図式化したもの。数字はサンプル数でプログラムと対応している。色の付いていないマスにはデータがないので値0を入れている。FilterのFFT処理はレイアウト上、矢印が引けなかったので、唐突に黄色く8サンプルで出てくるけど、上段のFilterの色がある3サンプルに値0の5サンプルを加えて8サンプルにし、FFTを通したもの。 Wikiにも原理の解説があったのでリンクしておきます。 http://ja.wikipedia.org/wiki/重畳加算法 入力信号のサンプル数の決定方法 下のプログラムではフィルタのサンプル数を元に計算している。 P=L+M-1 この式は、あるサンプル数の入力信号と、あるサンプル数のフィルタ信号を畳み込むと、必要なサンプル数が増えるので、そのサンプル数を求める式。 P: 出力の長さ(FFTに入れるサンプル数) L: FFTにかける入力サンプル数(全体の長さではない) M: フィルタのサンプル数 ここに分かっている数値を入れて、Lを導いてみる。 8=L+3-1 L=6 LはMよりも大きいという条件にしたが、自己流なので本来どのように計算するかは不明。いつものように少ない資料や式から、こじつけでプログラム化してます。今回は3サンプルなので、FFTは2^nなので、8サンプルもしくは、16サンプルが使えそうだ。ただ16サンプルだと無駄が多いので、この場合8サンプルが適当となる。ここまで決まったら、入力信号のサンプル数を計算する。 フィルタのサンプル数MとFFT用サンプル数Pが決まれば、上記計算式で入力信号は6サンプルで切ればよいことが分かる。FFTには8サンプル分入れて計算するので、残りの2サンプルには0を入れてしまう。同じようにフィルタもサンプル数

「DFTによる畳み込み」 記事と同じ内容をFFTでもやってみた。 内容は下図の通りで、1,1,0,0 の4つの入力信号と 0,1,0,0 のフィルタを、それぞれFFT(高速フーリエ変換)して周波数領域にしてしまう。そして周波数領域で複素数の乗算をして、iFFT(逆高速フーリエ変換)で戻すというもの。フィルタの内容は1サンプル遅れるディレイなので、結果は 0,1,1,0 となり、時間軸での畳み込みと同じ結果が出る。 FFTで重要なのはサンプル数が2^nであるということ。今回は4という少ないサンプル数で試している。 この記事の例はサンプル数が少なすぎて、現実感がない。実際の入力信号は、かなり長いのが普通だし、リアルタイムであれば、長さが不明ということになる。フィルタは普通数千サンプルぐらいで、長くても数万サンプル。この長さの違うサンプルを畳み込む必要がある。その方法は overlap-add method という畳み込み手法で、次回の記事で書いてみたいと思う。この方法を使えば、高速なFIRデジタルフィルタを実現できる。実際に多くのエフェクトで利用されている。Audacity の Equlaization や Convolution Reverb などがよい例。 FFT convolution 実験用ソースコード /* fft_test.c */ #include &ltstdio.h&gt #include &ltmath.h&gt /* FFT Cooley-Tukey型FFTアルゴリズム */ void fft(int *n, double *t, double *ar, double *ai) { int m, mh, j, k, irev, i=0; double wr, wi, xr, xi; for (j=1; j&lt*n-1; j++){ for(k=*n &gt&gt1; k&gt(i^=k); k&gt&gt=1){} if(j &lt i){ xr = *(ar+j); xi = *(ai+j); *(

まとめ作業は大事 プログラムをしないで半年経ってしまったが、昨年末から少しやり始めている。半年前にC言語を体系的に学習してBlogにアップしていたので、それを見ればちょっとしたプログラムはすぐに組めた。作ったプログラムは こちら 。後から見ることを前提に、まとめておくのは有効だと思えた。 C言語はやっぱり軽い 開発環境はMeadow(テキストエディタ)でコードを打って、 GCC を使ってコンパイルしてexeを作っている。このお手軽さに慣れてしまうと、もうJavaには戻れそうもない。 Javaの開発にはEclipseを使っていたが、あれが重い。開発環境はそれなりに巨大で動作がもっさりしている。それに起動も遅く、思いつきで何か書こうという気にならない。出来たアプリケーションのファイルサイズこそ小さいものの、起動は遅く、メモリは最低30MBぐらい使ってしまう。小さなアプリをたくさん作るという向きの言語ではないね。また最近ではJavaの脆弱性が頻繁に問題になっていて、Blog上にJavaAppletを載せるのもどうかという気がしてきた。いろいろJavaの使勝手の悪さが気になりだしてしまった。GUIも簡単にできるし、いろいろ便利で好きな言語ではあるのだが・・・ それに対してC言語 ANSI C89の範囲で作ったアプリは、すべてが軽い。特に最終的に出来たアプリはexeファイルサイズが小さい。GCC の「-s」オブションを兼用することで10KB前後（オプションなしだと50KB前後）になる。さらにネイティブなので高速に起動し、高速に処理をする。この差は大きい。やりたいことが結構重い処理である音声ファイルの加工だったりするので、C言語が適当だろう。今のところサクサク動くコンソールアプリ(CUI)で不満はない。そのうち wxWidgets などを使ってGUIアプリに手を出すかもしれないが、ちょっと試したところファイルサイズが大きすぎて魅力がない。それにGUIの必要性もあまり感じていない。GUIが欲しいアプリを作る場合は、やっぱりJavaを使うかな。 C言語で出来ないこと GUIは簡単ではない。Windowsであれば、結局Windowsのライブラリを利用するしかない。 wxWidgets を使っても同じこと。純粋なC言語の範囲ではコンソールアプリしか作

説明用としてDFTを使った畳み込みプログラムを作ってみた。本来はDFTを高速にしたアルゴリズムであるFFTで作るものだが、流れを確認するために簡単なDFTにしてみた。 式とかはこちらを参考にした。 http://ja.wikipedia.org/wiki/畳み込み 下のプログラムでは、入力サンプルは4つで{1,1,0,0}というもの。フィルタも4つで{0,1,0,0}というもの。これは1サンプル遅れて出力されることを意味する。つまりディレイ。出力結果は{0,1,1,0}となる。この例は単純であるが、実際には周波数成分のコントロールや残響のコントロールなど幅広い応用が考えられる。 これを実現する方法で一番シンプルな方法は時間領域で畳み込むこと。しかし、データ量が増えると処理速度の問題が出てしまう。そこで、下図のように一度DFTで周波数領域に置き換えて、そこで演算を行い、その結果をiDFTで戻して出力をするという一見面倒な処理をするのだが、FFTのおかげで時間領域の畳み込みよりも高速に処理することができる。ここでは学習向けにDFTを使っているが、本来はFFT。FFTによる畳み込みはFIRの各種デジタルフィルター、EQ、HPF、Convolution Reverbなどで広く利用されている。 同じ内容をFFT でもやってみた。 DFT convolution 実験用ソースコード このプログラムは数式をそのままにして自己流で作ったもの。参考にしているプログラムはないです。ということで間違っている可能性あり。また上から順に説明するために処理を関数にしてないです。複素数計算のところは怪しい。本来はもっとスマートに書くのだろう。またcomplex.hの利用も考えたけど、直接いじっている感じがしなかったのでやめた。またやたらと配列を使っているが、これも実用的なプログラムの場合はポインタを使う。まぁいろいろ突っ込みどころの多そうなサンプルです。そのかわり、とても素直で理解しやすいと思う。 /* dftcon.c */ #include &ltstdio.h&gt #include &ltmath.h&gt /************************************************

2012年1月の記事にプログラム追加。Javaで作ったプログラムをC言語に移植し、iDFTで信号を戻す部分を追加した。 DFT（discrete Fourier transformation）とは離散フーリエ変換といって、離散信号の時間領域の波形を周波数領域に変換する手法のこと。 下の例では波形をDFTによって分析してみる。 上記波形をDFTで周波数領域に変換すると、下のように1000Hzと10000Hzの音が含まれていることがわかる。 DFTは音関連のプログラムを作る上で、とても基礎的なことなので、プログラムを作りながら確認してみることにした。 基本となる式はこれで、H(k)は最終的に欲しい各周波数の振幅。h(n)は分析する各サンプル。Nはサンプルの数。この逆の操作となる逆フーリエ変換は下のようになる。 上記の式のままでは使えないので、使いやすいように下のように式を変形する。下付き添字のrは実部を表し、iは虚部を表す。数式はあまり詳しくないので、書き方が間違っているかもしれないが、動作としては一番下のプログラムソースを参考にしてください。 逆変換は下のようになる。 フーリエ変換のプロセスは結構面白いので、そのうち子ども向けに解説しようかと思っているが、今回は流れだけ簡単に。 まず、周期性のある波形は、その整数倍のsin、cos波に分解できるという、フーリエの理論がある。具体的には、調べたい波形に整数倍の sin、cos波を掛けて積分すると、各周波数成分が、どれぐらい含まれているかが出てくるという仕掛け。含まれていなければ0で、含まれていれば、振幅のレベルが出てくる。このプロセス自体はシンプルなのだが膨大な計算量なのが欠点。より高速に処理するためのFFT(高速フーリエ変換)という手法もあるのだが、今回はシンプルなDFTをまず作ることにした。上記の式をプログラムで実行するのは簡単で、入れ子のfor文数行で完成してしまう。ただ式に虚数 i （電気ではj）があるので複素数を扱う。周波数スペクトルに使う場合は、実部と虚部の数値をそれぞれ２乗して足した平方根が、周波数スペクトルの振幅スペクトルとして利用できるようになる。式はこんなかんじ。 このXnはリニア表示なので、さらにAudacity

57-63日目 初録音 曲を最初から最後まで弾けるようになったので録音してみた。パソコンに手作りマイクを挿してAudacityで録音するという簡易的な録音。ギターを2本重ねたいと言うので、なぜかダブリングを行う。これは2本同じ演奏をして音を重ねてしまうという方法で、ぴったり演奏すると2本のギターが鳴っているというよりも、音に厚みが出る効果がある。録音後聴いてみると意外とそれなりに合わせていた。はじめっから、なかなかやるなぁ。 歌の録音も ギター弾きながら、歌えないというので、ギターを録音してから、それを聴きながら歌ってみるという方法で録音。録音が終わって始めて自分の歌を聴いてみる・・・ 相当ショックのようでした。「自分の思っている声とは違う」もしくは「自分はもっとうまく歌っているつもり」のどちらかだろう。客観的に判断するのに録音はよいのだが、はじめは抵抗あるよね。プロでも自分の声を嫌だと言う人多いし。そんでもって、もう歌わないようなことを言っていたけど、後日普通に歌っていたので、立ち直ったみたい。 8週目＜     ＞10週目以降 ギター練習

先月LADSPAプラグイン 同じようなもの を作ったのだけど、より手軽にNyquistでも作ってみた。Nyquistはインタプリタなので、初期設定を簡単にいじれるよさがある。8つの任意の周波数とゲインのサイン波を合成して生成するプラグイン。AudacityではGenerateタイプになる。この種類のプラグインは、はじめて作った。 インストール 使用するにはテキストエディタに下記ソースコード内容をコピー。名前はnamagi_sinewave.nyとして保存。Audacity/Plug-Insフォルダの中に入れれば使えるようになる。 使いかた AudacityのGenerateメニューからNamagi: SineWaveを選択すると以下のウィンドウが表示されるので、8つの周波数とゲインを調整し、生成する時間を入れてOKを押せば、合成波形が作られる。 上記の設定のまま実行すると、以下のような複雑な波形が10秒作られる。 ソースコードはとても単純。使っているのはNyquistに用意されたものだけなので、あまりプログラムとは呼べないレベル。stretchで生成する時間を指定して、あとはひたすらsineを足していくという内容。 Namagi SineWave Nyquist ソースコード ;nyquist plug-in ;version 1 ;type generate ;name "Namagi: SineWave ver.130201..." ;action "Namagi: SineWave ver.130201..." ;info "" ;control time "length(sec)" real "" 10 0 600 ;control freq1 "freq1" real "" 100 0 22050 ;control gain1 "gain1" real "" 0.2 0 1.0 ;control freq2 "freq2" real "" 200 30 22050 ;contro