FFT - みる会図書館


検索対象: 月刊 C MAGAZINE 2002年3月号
12件見つかりました。

1. 月刊 C MAGAZINE 2002年3月号

オーディオの周波数解析手段として一般 に広く用いられている短時間フーリエ変換 ですが , 実際に存在しているスペクトルを 厳密に描写しているわけではないというこ とは覚えておいてもいいと思います。だい たいの傾向が一致するという , あくまでも 近似処理なのです。しかしこれはたとえば フーリエ変換が間違っているなどといって いるのではないことに注意してください。 短時間フーリエ変換の結果をその区間の物 理的なスペクトルだと考えると少々無理が 出てくる , というだけのことです。 DFT と FFT の導入 では今度は , 実際に周波数解析を実装し ていくまでを見ていきます。 式⑩で与えられている短時間フーリエ変 換は連続関数 ( アナログ信号だと考えられ る ) に対してのフーリエ変換であるため , このままではデジタル信号には適用できま せん。そこで短時間フーリエ変換の離散デ ータ用パージョンである DFT(Discrete Fo urier Analysis) を使うのが普通です。この D を式 2 に示します。 先ほどの時刻 tl ~ あという範囲が , 単にデ ータ数Ⅳの範囲として変更されています。 しかし , 実はこのの式が問題なのです。 この式をそのまま計算するとかなり重い処 理なのです。そこで PC 上などで周波数解 析を行う際には , FFI 、 (Fast Fourier Analys (s) と呼ばれるアルゴリズムを使うのが普 通です。 F はバタフライ演算という演算 方法によって D 矼をアルゴリズム的に最適 Fig. 26 FFT 演算前後のデータフロー x はト Ex 回ビ 化したもので , 速度は劇的に改善されてい ます。だいたいの感じとしては , D ドでは その結果を得るのに PC が無反応になり , し ばらくしてやっと結果が表示されていた ( 笑 ) ものが , F ド I 、ではポンポンリアルタイムに 追従できるくらいのレベルになります。 FFT のアルゴリズムの内部などに関して は , 参考文献そのほかを参照してください。 FFT を使った信号解析 さて , それでは F を使って実際に時間 軸上のデータを周波数軸上のデータに変換 するまでを解説します。この手順は以下に 示すとおりです。 ①サンプルの信号の切り出し ②切り出した信号に窓関数を乗じる ③ FFT に実数部テータとして代入 ( 虚数 部はすべてゼロとする ) ④ FFT 演算 ⑤結果を適宜パワースペクトルに修正 この手順を直感的につかめるように示し たのが Fig. 26 です。 F 、は通常は複素演算 であるため , 図のように実数部と虚数部に 分けて同じ長さ N のデータを入出力するの が普通です。実数部に実信号データを入れ , 虚数部はすべてゼロとします。ただし F ロ、 の結果としては虚数部もゼロでなく値を持 つことに注意してください。また出力結果 は , 実数虚数とも左右対称です。 これは F ドがゼロ ~ サンプリング周波数 までの周波数解析ができるため , サンプリ ング周波数の半分のナイキスト周波数を境 にイメージ信号が現れるために起こる現象 実数部 ( 窓関数を乗じたⅣサンプルの信号 ) 0 p 9 0 0 p ′ 0 です。 切り出す入力データ数 Nt こは約束があり ます。 F ではこの Nt こ任意の数は許され ず , 256 , 512 といった , 2 " の数でなくては いけないことになっています。 F ド I 、演算内 部ではこのⅣが 2 。であることを利用して効 率的な計算をしているためです。また , F 後のスペクトルの周波数間隔 4 / ・は , サ ンプリング周波数までに N 組の実数部と虚 数部のデータが出力されることから , サン プリング周波数を Fs とすれば , 式で与え られます。こで与えられる」、 / ・が周波数 の最小解像度になるため , 目的に合うよう に N を決めていくことになります。リアル タイム FFT の場合 , は 512 ~ 2048 くらい が相場でしよう。 また , F ドの出力結果のうち実数部を [k] , 虚数部を l[k] とすると , そのパワー スペクトル P [ k ] は式で与えられます。 つまり , R[k] とⅡ k ] を k 番目の複素数の 実数部と虚数部と考えた場合 , パワースペ クトルはその複素数の大きさに相当します。 この PCk] をプロットすると Fig. 1 のような スペクトル情報になります。 FFT による解析部分の実装 周波数解析の話の最後に , F ドを使った 周波数解析のソース部分を見ていくことに しましよう。 List 3 を見てください。 FFP の 計算はこの CFilterStudioDlg::CaIcFFT 関数 の中で行っています。まず , CStreamMan ager::GetFITData 関数を使って , 現在演奏 している付近のデータを FFI 、 SIZE サンプ ルだけ左右チャンネルそれぞれに抜き出し ます。これに備えて CStreamManager は「最 近の」データをキャッシュしています。 の後あらかじめ計算してその値を m ー pWin dow に格納しておいた窓関数をデータに乗 0 0 0 、 0 実数部出力 b00900d 9 0 0- 2 、〇一 0 ◇づ 虚数部出力 ーノ 2 ド FFT 演算 28 虚数部 ( すべてゼロ ) C MAGAZINE 2002 3 0-0 -0-0 一〇一 0 -0 -0-0 ー 0 ー 0 -0 ・ 0

2. 月刊 C MAGAZINE 2002年3月号

LIST- 3 F FT の計算 void CFi lterStudioD 地 : : Cal cFFT ( ) / / 実数部にはデータ fo て ( 土 = i く FFT_SIZE; i 十十 ) { / / 窓関数の乗算 int 土ー m—Lock. Unlock( m—StreamManager. GetFFTData ( m—pLeftData , m—pRightData , FFT—S I ZE m—Lock. Lock ( / / F のためにデータの取得 fo て (i = の土く FFT_SIZE / 土十十 ) { / / 実数部と虚数部からパワースペクトルを計算 m—FFT. Comp ー exFFT ( m_pRightFFTData m—FFT. ComplexFFT(m—pLeftFFTData); / / 複素 F 計算 m—pRightFFTData[i * 2 十 1 ] = 0.0 新 〃虚数部はゼロ m—pRightFFTData[i * 幻 = m—pRightData[i] * m—pWindow[i]; 〃実数部にはデータ m pLeftFFTData[i * 2 十 1 ] = 0.0 / / 虚数部はゼロ m-pLeftFFTData[i * 2 ] = m-pLeftData[i] * m-pwindow[il; m—pRightFFTData[i * 2 十 1 ] * m—pRightFFTData[i * 2 ] 十 m-PRightData[i] = sqrt(m-pRightFFTData[i * 2 ] * m—pLeftFFTData[i * 2 + 1 ] 潺 m—pLeftFFTData[i * 2 十 1 ] * m—pLeftFFTData[i * 幻十 m—pLeftData[i] = sqrt(m—pLeftFFTData[i * 2 ] * m-pRightFFTData[i * 2 十 1 ] / / 半分のみを描画データとして渡す / / F の結果は FFT-SIZE / 2 で左右対称なので、 m—wndScreen. DrawFFTResult(m—pLeftData, m—pRightData, FFT_SIZE / 00 Fig. 27 周波数軸上での操作 時間軸データ 元の 加工された 周波数軸データ シフトさせる 周波数軸上で 処理後 処理前 Fig. 29 時間軸データ テフり FFT 逆 F 阿 テータ操作 何らかの 2 算します。その後 , CFFT::ComplexFFI 、関 数で実際の F ド r 言 t 算を行うわけですが , こ の関数は入出力の実数部と虚数部がそれぞ れの独立した配列に分かれておらず , 1 つ の配列に交互に「実虚実虚」と並ぶ仕様に なっていることに注意してください。 FFI 、 の計算が終わったら , 式@で示されるパワ ースペクトルの計算をし , F ド I 、の表示クラ スにデータを渡しています。 ところで , バタフライ演算に代表される F 、の内部計算はこれまたたいへんなノウ ハウが必要な部分で , 実装によって速度は 大きく変わります。今回採用した CF ドと いうクラスは , 大島氏 ( 参考文献 U Ⅲ滲照 ) のソースを筆者が C + + で薄くラップしたも のです。 F ドの内部計算に興味がある方は 別途文献などを参照してください。 なお , F ド I 、の更新速度は FFT_INTERVAL で指定しており , デフォルトは 30 ミリ秒で , Fig. 28 ピッチスケーリング 1 秒間に 30 回程度 F の計算と表示をする ことになります。同時に走らせる機能にも よりますが , 最近の CPU ならば問題なく動 作するはずです。どうしてもきつい場合は FFP_INTERVAL を 100(100 ミリ秒 ) という 特集 1 デジタルフィルタソフトウェアによるオーディオデジタル信号処理のレシピ 29 ピッチスケーリング前後の波形と周波数成分 時間は変わらず周波数成分のみ変化 ように値を大きくしてみてください。 このような音楽再生をリアルタイムに追 従するような周波数解析は , 眺めているだ けでもけっこうおもしろいものです。最近 のオーディオアプリケーションは , デザイ ンを凝らしたうえでイルミネーション代わ りに表示させているものもあります。たい へんなぜいたくができるようになりました。 第を T を使った ー。 . フェクト処理 こまでで , F ロ、を使った周波数解析機 能の実装について解説しました。今度は F ド r を単なる解析のツールとしてではなく , もっと積極的にエフェクトとして活用して いくことを考えていきたいと思います。 ピッチスケーリング 実は F ド I 、のべースになっているフーリエ

3. 月刊 C MAGAZINE 2002年3月号

うなシステムをフェーズボコーダといいま す ( アナログ時代からあるチャンネルボコ ーダとは意味が少々異なる ) 。フェーズボ コーダを使うと各正弦波の振幅情報や位相 情報が独立に操作できるため , ピッチスケ ーリングに代表されるようにたいへん自由 度の高い処理を実現することができます。 ピッチスケーリングの実装 ピッチスケーリングの構成がわかったと ころで , この機能を実際に実現するソース の解説に移ります。 List4 を見てください。 こまで解説してきたピッチスケーリング は正逆の FFT 計算以外はすべてこの C PitchScaIeUnit クラスで行っています。なお , LR のチャンネルごとに , このクラスのイン スタンスが必要です。 CPitchScaleUnit::Operate 関数がピッチス ケーリングのメイン関数です。パラメータ に入力サンプル数 , 入出力用パッフアを指 定します。 F ド / 、、渡すデータサイズの制約 から , スケーリング処理は m nFFTSize ( デフォルト 2048 サンプル ) ごとにしか行う ことができません。このため毎回の入力デ ータサンプル数がいかなる値であっても , 内部に 2048 サンプルたまるまではデータ処 理を待つようなキャッシュ処理をしていま す。このキャッシュ領域となるのが , m_p InputBuffer です。 さて , キャッシュサイズが F ドできるだ けの大きさになったら , まずはコサイン窓 を掛けてデータをキャッシュから切り出し ます。ここでのコサイン窓は前述した窓関 数の一種ですが , サイドロープ低減よりも むしろできるだけスムーズに両端を絞り込 んで , 不連続点が出ないようにするために 使っています。 m_pFFTBuffer が FFI 、に渡さ れるデータバッフアですが , この m_pFFT Bu 幵 er にはやはり「実虚実虚」と実数と虚数 が交互に並びます。 m-pFFTBuffer の準備 ができたら , 正 F ド r 言 t 算用オプジェクト m FFTForward に渡し FFT を行います。この 結果は真の周波数推定のため , GetTrueSp ectra 関数に渡されます。次に ShiftTrueSpe ctra 関数で , 得られた真の周波数によるス 32 C MAGAZINE 2 側 2 3 ペクトル列を所定の量だけスケーリングし ます。シフトされたスペクトル列は , Crea teShiftedVector 関数で再び実数部と虚数部 に分けられ , メンバ変数 m_pFFTBuffer 内 に準備されます。そして最後に逆 F ド言 t 算 用オプジェクト m FFTReverse で時間領域 に戻されるようになっています。 時間軸データに戻されたデータは , 出力 キャッシュ m—pOutputBuffer に再びコサイ ン窓を掛けて足し込まれます。なお , Fig. 32 の近接した F ロ、フレーム間の距離は m ー n OffsetSize にあたり , デフォルトでは 512 サ ンプル ( 2048 / 4 ) となっています。なお , 入 出力キャッシュではフレーム間で Fig. 34 の ような重ね合わせが行われていることに注 意してください。この重ね合わせのために ゲイン ( 振幅方向 ) も補正しています。 Fig. 34 では表現しきれませんでしたが , 1 つの フレーム操作が終わったら入出力キャッシ ュ内のデータを memmove 関数を使って m nO 爪 etSize ぶんだけ左方向にシフトしてい ることに注意してください。こうしないと フレーム処理が進むにつれてどんどん入出 カキャッシュの配列のインデックスが大き くなっていくことになり , 固定長キャッシ ュではまかなえなくなります。 個々の関数の中もざっと見てみましよう。 GetTrueSpectra 関数内部では , まず実数部 と虚数部データから , 1 本 1 本に分解された 正弦波の振幅と位相を計算します。この位 相情報を元に , 前回のフレームとの位相差 を計算し fDeItaPhase に代入します。次にこ の位相差からフレーム間の時間差によるず れのぶんを引き , F ドグリッドと真の周波 数との差による位相ずれ fOffsetPhase を求 めます。なおここでの位相ずれはマ ~ + の間でないと , 隣接するの周波数グリッド を超えてさらに隣のグリッドにまで補正が 影響してしまいます。このため正規化処理 を行い , 位相ずれを - ~ + の範囲に正規 化します。この処理は一見せつかく得られ た位相差を捨ててしまっているように見え ますが , 正弦波はもともと - ~ + で周期 的な波形なので , これでかまいません。そ してこの位相差を元に , 真のスペクトルが 幅 1 に正規化しているグリッドからどれだ けずれた位置にあるかを fO Ⅱ set に格納しま す。最後に fOffset とグリッドインデックス i を足したものにグリッド間の周波数差を 乗じて , 真の周波数 fTrueFreq が求められ ます。 ShiftTrueSpectra 関数ではスペクトルの 周波数方向のシフトを行います。振幅情報 はそのままにして , 周波数情報のみを fSca le だけシフトして , メンバ変数に格納して います。 CreateShiftedVector 関数は振幅と真の周 波数に分かれたスペクトル情報を , 逆 FFT 用に再び実数部と虚数部に分配する関数で す。真の周波数 (m_pSynthesizedFreq [i] ) と最寄りの F ロ、の周波数 (i * fDeltaFreq) とのずれのぶんだけわざと位相をずらして います。フレーム間の時間差による位相差 も加えておきます。 外部からのパラメータ設定について解説 します。この CPitchScaleUnit クラスを使う 側から , 初期化時に各種のパラメータを設 定する関数が CPi 忙 h eU ⅲ t : 面 it 関数です。 パラメータに FFT サイズ , サンプリング レート , オーバラップ数を設定します。 F ロ、サイズは 2 。でなくてはならず , 1024 か 2 048 が妥当です。音質の観点からはできる だけ大きくしたいのですが , そのぶん重く なります。サンプリングレートにはサンプ リング周波数を Hz の単位で指定します。オ ーバラップ数は , Fig. 34 に示すように , FF T のフレームをどれだけ隣と重ねるかを指 定します。 Fig. 34 の場合 , フレームの 1 / 4 が重なっているので 4 を指定します。これ また音質のためにはできるだけ大きい数を 指定したいのですが , やはり重くなります。 なおオーバラップ数は F ドサイズを割り切 れる数にします。 FFT サイズは 2 。なので , オーバラップ数も 2 , 4 , 8 , 16 といった 2 。で ある必要があります。ただし内部計算での 多少の誤差を容認できるなら , オーバラッ プ数は 2 。でなくてもいちおう大丈夫で , 音 もそんなに違和感はありませんでした ( た だしオーバラップ数が F ドサイズを割り切 れない数になるのは , 理論的には少々おか

4. 月刊 C MAGAZINE 2002年3月号

処理の重さに直結しますし , あまりに幅が 広いとせつかく分割したフレームが原音の とおりにくつついてしまう ( 戻ってしまう ) こともあるので , 注意が必要です。デフォ ルトでも処理がきつい場合は , この値を小 さくしてみてください。 式⑩の符号評価は Sign 関数で行っていま す。この関数は実質上符号ビットを見るだ けでよいので , アセンプラを使うなどして いろいろ試してみました。しかしリストに 示すようなゼロとの単純比較が意外に速く 見通しもよいので , 結局この方式を用いて います。 最良点 k, 。が見つかったら , 今度はフレー ムの連結を行います。フレームの連結はソ ースに示すようにフレーム同士がオーバラ ップする領域と , そうでない領域に分けて 行います。もちろん , オーバラップする領 域ではフレーム間にクロスフェーダ処理を 行います。フレームの連結が済んだら , nS a ぶんだけ入力キャッシュを左にずらして おきます。こうすることにより , 入力キャ ッシュの先頭が常に次のフレームの先頭に なります。 ISM 処理が済んだら , 出力バッフアに出 カサンプルをコピーします。出力キャッシ ュである m_pLeftY と m_pRightY から memc py 関数を用いてコピーします。なお , スケ ールレートの変更は CTimeScaleUnit::SetR ate 関数で外部から行います。 Fig. 35 に示すように , サンプルアプリケ ーションでは公開パラメータとしてスケー ルレートを 0.5 ~ 2.0 くらいの範囲で連続可 変できるようにしています。 今回はこれらのどのレートでもそれなり に聞こえるようにという基準で , S , S , N といったフレームパラメータを試行錯誤し てみました。これらは最終的には順に nSa = 2000 / fRate, nSs = 2000 , nN = 2800 という値 となっています。これらの値は 44.1KHz 専 用で , ほかのサンプリングレートの場合は 再調整が必要でしよう。なお , nSs は毎回 固定で , nSa のほうをスケールレートに応 じて可変にしていることに注意してくださ い。このほうがよい結果になるようです。 40 C MAGAZINE 2002 3 しかし fRate がどこかの値に固定の場合 は , 個々のケースに応じてもう少し音質的 に好ましいパラメータが見つかるかもしれ ません。このあたりはソース依存だけでな く個人差もあり , 一概にどういうパラメー タ値がよいということは難しいでしよう。 もしもこのままで不満が残るようなら , カ ット & トライでよいパラメータを探してみ てください ( この爲 M を試す場合にはリリ ースビルドで行ってください。デノヾッグビ ルドだと EM - TSM の評価計算に時間がかか るようで , ブップッと音が切れる場合があ ります ) 。 PART3 のまとめ 本章では爲 M と呼ばれる , 速度変更アル ゴリズムの理論から実装までを見てきまし た。今回取り上げた 0 はをベースにした操 作は音声合成のための要素技術として発展 してきた感があり , どちらかというと単一 音声向けです。そのためたくさんの楽器が 入っているような今日の超 Hi-Fi オーディオ の基準で考えると , クオリティ的には多少 もの足りない気がするかもしれません。し かし筆者が初めて SOLA を試したときには 「こんな処理でもテンポを変更したように 聞こえるんだ ! 」という部分でおおいに驚 いた記憶があります。みなさんもぜひサン プルアプリケーションで試してみてくださ い。かなり強引な処理をしているわりには 「意外に」音がよいことがわかっていただけ るでしよう。スケールレートが 1 から離れ るにつれてどんどんクオリティはきつくな ってしまいますが , 逆に士 20 % 程度なら音 質はかなり満足できるレベルで , 複雑な波 形のオーディオソースにかけてもほとんど 違和感はありません。 また何かの会話を録音したものを後で 倍速で速聴したりするときなどには , それ ほどのクオリティはいらないでしようから , 逆にこのエフェクトはかなりの威力を発揮 しそうです。まるで早ロ言葉のように聞こ えることでしよう。 さてこれで , ・時間軸のみの操作 ( PART3 ) ・周波数軸のみの操作 ( PART2 ) ・時間軸と周波数軸同時の操作 ( 2002 年 の 6 F 阿と Wave let 本稿で FFT をトピックの 1 っとして取り 上げていますが , 実は FFT はけっこうクセ モノです。 FFT は , 変換結果が間隔ム f = F / Ⅳでリニアスケール上に等間隔に並ぶ変 換です。ところが , 実際の信号処理で便利 なのは , リニアスケール上ではなく , ログ スケール上で等間隔に並ぶような変換なの です。これは PART2 で , ピッチスケーリ ングのスケーリング操作が , ログスケール で等間隔になるように掛け算シフトをした ことでわかっていただけるでしよう。フェ ーズボコーダで各成分に分けたような操作 を各バンドバスフィルタでの分解と考える し FFT は各フィルタが一定の Q ではなく 高域になるにつれて Q が上がってくる変換 といえます。 変換結果のスペクトルがログスケール上 で等間隔に並んで Q が一定なのが , 本誌で もたまに取り上げられる wave 回変換です。 特筆すべきは , Wave 回変換の解析結果は , 個々に独立した要素となっているという性 質です。もしも FFT で変換した 1 つのスペ クトルを勝手に消去してしまい再度逆 FFT で戻せば , たいていの場合区間の両端で不 連続点が発生するという困った問題が起こ ります。これは FFT では解析結果の正弦波 が区間内 ( 実は区間の外でも ) で一様に連 続していると暗黙のうちに仮定しているた めに起こる現象です。これに対して Wavel et 変換は , 周波数と場所とでもいうような 2 つのパラメータを使うことで , 区間外は おろか区間内でも各成分が局在化している という前提のもとで解析を行います。この ため , ある成分を勝手に抜いてしまっても , 不連続点のような不整合が起こりにくくな っています。 Wavelet が注目され始めたのはわりと近 年で、現在もまだ研究が盛んに行われてい る段階です。今回紹介したピッチスケーリ ングの Wavelet 版のようなものもすでに存 在し , これからが楽しみな技術といえるで しよう。

5. 月刊 C MAGAZINE 2002年3月号

P A R T FFT を使った信号処理 本章では , FFT ( フーリエ変換 ) を使った信号処理について 見ていきます。前半では FFT を使った周波数解析の理論に ついておさらいしながら , 最終的には WAVE ファイルの スペクトルアナライザを実装します。そして後半では FFT をもっと積極的に信号処理に応用する例として , ピッチス ケーリングについて解説します。ピッチスケーリングは再 生時間を変えずに , そのピッチのみを変更するテクニック です。音楽でいう移調のようなものと思えばよいでしよう。 こうした操作を行うためには信号の周波数軸上での操作が 基本になりますが , 時間軸上の信号データを周波数軸上の データに変換できるのがまさにフーリエ変換なのです。 FFT による周波数解析 周波数解析 ( スペクトル解析 ) とは , どん な周波数の信号がどれくらいのレベルで含 まれているかを解析する手段です。最近は MP3 プレイヤでもこの機能を持つものが多 くなり , 非常にポピュラーな解析になって います。本章ではまずこの周波数解析機能 の実装について理論から解説します。 フーリエ変換 Fig. 1 のサンプルアプリケーションは , リアルタイムで周波数解析を行う機能がつ いており , その結果を下側のグラフのよう なエリアに表示することができます。べー スが鳴れば低い周波数成分のレベルが上が り , ストリングスが鳴れば高い周波数成分 のレベルが上がります。水平方向右側にい くほど高い周波数になり , 垂直方向上側に いくほど , レベルの高い周波数成分となっ ています。 しかしよく考えてみてください。通常の WAVE ファイルや MP3 などのオーディオフ ァイルは , PARTI で示したように時間軸上 のデータです。これに対して , Fig. 1 のグ ラフは周波数軸上のデータ分布です。した 26 C MAGAZINE 2002 3 がって , Fig. 1 のようなスペクトルデータ を得るには , もともと時間軸上のデータと して提供されている信号を周波数軸上のデ ータに変換する手段が必要なのです。この 変換を提供するのが , 有名なフーリエ変換 です。フーリエ変換は文献も多く , また本 誌でもしばしば取り上げられているトピッ ここではポイントを絞って軽く クなので , おさらいしていきましよう。フーリエ変換 自体は純粋に数学的な操作で , 式⑩のよう に表現されます。 フーリエ変換はさまざまな分野で見かけ ますが , もしも信号処理の世界で解釈する なら , 先にも触れたように「時間軸上のデ ータ x ⑦を角周波数軸上のデータ F ( の ) に 変換する処理」ということになります。 れは式⑩が実質的には元のべクトル x(t) を各周波数成分 ( スペクトル ) に分解する こと ( 各成分が直交しているため , べクト ルの成分分解にあたる ) に対応しているた め可能になる機能です。 しかしここで注意すべきは , 式⑩は時間 軸での無限積分となっている点です。時間 軸での無限積分ということは , 実際の物理 に則して考えるなら , さしずめ「宇宙の始 まりから終わりまで」とでも表現すべきも のであって , 数式では気楽に表現されてい ても実際はなかなか一筋縄ではいきませ ん。 そこまで極端に考えずに , 仮にこの積分 範囲が音楽 1 曲の始まりから終わりまで ( それ以外はゼロと考える ) だとしても , これまた問題です。 1 曲全体を解析してや っとひとまとめの周波数成分に分解できた としても , あまり実用性がありません。オ ーディオアプリケーションでは音楽の特定 部分だけを局所的に解析するのが普通だか らです。そしてこの解析場所を音楽の再生 に同期させて進めていき , その瞬間ごとに 音楽に含まれている周波数成分を観察する わけです。 しかし , フーリエ変換を普通にそのまま 考えるだけでは , 積分区間の問題のために 1 曲全体に対して 1 セットの変換結果しか 得られず , これは「曲全体を通してどの周 波数がどれだけ出ているか」という静的な 解析になってしまいます。刻一刻と変化す る音楽信号に含まれる周波数に応じてビョ コピョコとグラフが動く Fig. 1 のような動 的な解析は , そのままのフーリエ変換では できないのです。 短時間フーリエ変換 このような問題に対して考えられたのが 短時間フーリエ変換と呼ばれるもので , 式 ⑩に示される変換です。 式⑩の積分範囲が乙から t2 に変更された だけですが , その意味は大きく変わります。 先ほどは無限大の区間必要だった信号デー タが , 今度は時刻 ~ しの範囲でのみ得られ ればよいからです。つまり積分に必要な時 間区間が短くなったことに相当するため , 短時間フーリエ変換と呼ばれるわけです。 しかしもともとは無限積分だったものを 勝手にある範囲に丸め込めてしまうため , やはり不都合は起こります。 Fig. 22 を見て ください。実信号が図の上段だとして , ⑩ F(O) = / ズいレーノ " 市 ⑩ F(O) = / ズはレーノ " 市

6. 月刊 C MAGAZINE 2002年3月号

変換には逆変換があります。 F ロ、でも当然 逆変換が存在し , 逆 F Ⅵ、と呼ばれています。 逆 F ロ、は正 F の文字どおり逆で , 周波数 軸上のデータを時間軸上のデータに戻すこ とができる変換です。 カンのよい方なら , この性質を使うとい ろいろおもしろいことができることに気づ くかもしれません。 Fig. 27 を見てください。 元のデータを F ドでいったん周波数軸上の データに変換し , このデータに対して何ら かのデータ操作をします。その後に逆 F ド で再び時間軸上のデータに戻せば , 時間軸 上での操作では扱いにくいようなさまざま なエフェクトを実現することができます。 今回はこうしたエフェクトの一例として , ピッチスケーリングを取り上げます。ピッ チスケーリングは Fig. 28 に示すような操作 で , 周波数軸上で各スペクトルを上下にシ フトさせる処理です。この処理を行うと , Fig. 30 1 組の実数と虚数の組で正弦波を表 虚数軸 ( 日図川 k]) の 日図 2 + / 図 R[kl 実数軸 元の波形の周波数成分 ( ピッチ ) を上下に シフトするエフェクトができることになり ます。特筆すべきは Fig. 29 に示すように この処理の前後では音の時間的な長さは変 わらないことです。 音楽で考えるならば , これは元データを 移調したことにあたり , 事実この処理を実 装すると , たとえば WAVE ファイルを 3 度 上げて再生したり , 1 オクタープ上げて再 生したりといったことが簡単に実現できま す。原調を移調してやれば , カラオケ練習 アプリケーションができるかもしれません。 こうした処理は時間軸上でも不可能ではあ りませんが一般に困難で , まさに周波数軸 上ならではといえるでしよう。 正弦波への分解 こで F ドで得られる解析結果を , もう 少し詳しく見ていきましよう。フーリエ変 換はもともとは与えられた信号をさまざま な周波数の正弦波に分解する変換です。 F FT でもこの性質はそのまま引き継がれて おり , 解析結果として得られる実数部と虚 数部のデータがワンペアになって , 1 つの 正弦波を記述する形になっています。 この様子を Fig. 30 に示しました。 k 番目 の実数と虚数のペア (R[k],I[k]) で複素平 面上に点をとり , この点の原点からの距離 P[k]= 心 k ] 2 + Ⅱ k ] 2 ( パワースペクトルに 相当 ) が記述される正弦波の振幅成分を示 します。また , 図中に / で示した ' のペア のアングル , すなわち " 。ね n ( - 々 [ 訂 ) で Fig. 31 実際のスペクトル位置と FFT の解析結果のずれ 最寄りのグリッド上に出てきてしまう FFT の解析では本当の位置ではなく 本当のスペクトル位置 Fig. 32 位相差計測法 F FT その正弦波の位相成分を記述します。つま り , 数式で書くなら , 式になる正弦波が 描写されていることになります。早い話が 式は直交座標から極座標に変換している のです。結局 , サンプル数 N の F ド I 、では , N 本の上記正弦波に信号が分解されているこ とになります。 ところが問題は式の、 / 、 k です。このム ( 周波数に相当 ) は間隔 4 / 、で並んでいる周 波数グリッド上の k 番目の値なので , 式 になります。しかしこの周波数グリッドは F の解析によって自動的に決まる間隔で あって , 実際に信号の中に存在しているス ペクトルの周波数は必ずしもこの値とは一 致しません。 Fig. 31 を見てください。 つまり , ピッチスケーリングのようなシ ビアな周波数情報が必要な場合には「 F ド I 、 の機械的な解析結果」をそのままうのみに して処理するだけでは周波数精度が足りず , 「真の周波数によるスペクトル」に補正した うえで Fig. 28 のような操作を行わなくては ならないのです。 この意味からは , 本章前半で解説してい る周波数解析は信用ならないといえるかも しれません。しかしオーディオ帯域全域に わたって周波数分布の「傾向」を見るよう なラフな解析であれば本章前半の方法でそ ④は ] 2 + / は ] 2 sin 2 な + arctan こく近い場所で 2 回 FFT を行う F FT 「位相の差」に着目する 同じた番目のスペクトルの 30 C MAGAZINE 2 2 3

7. 月刊 C MAGAZINE 2002年3月号

プログラミング技術情報誌・ C マガジン C 0 N T E N T S デジタルフィルタ ソフトウェアによる オーディオデジタル信号処理のレシピ PARTI デジタルフィルタ PART2 FFT を使った信号処理 PART3 特殊な信号処理 C/C + + /Java 実力チェック プログラミング期末試験 C 言語編・・・出題・解説 : きだあきら C + + 編・・・・・出題・解説 : 大城正典 Java 編・・・・・出題・解説 : きだあきら GetInto C WorId ! ! ー c 言語入門講座 く最終回 > 入門レベル補足小薗三典 / 中井信会 C プログラマのための C + + 入門 実践 C + + ゼミナール く最終回 > C + + アプリケーション作成のポイント吉野興 スタートアップ Java—Java 言語事始く最終回 > コレクシ = ン毛呂宗夫 Enjoy PerI Programming モジュールを活用しよう く第 14 回 >YukiWikiMini/Template の改良結城浩 なあ Ruby を読もうじゃないか = ー = 。。 , 。ー = 。。 画像処理を極める ア丿レゴ、リス、ムラボく第 30 回 > 画像圧縮アルゴリズム⑨動的辞書を利用した圧縮昌達 K'z 特集 1 小川要 特集 2 13 42 86 93 98 107 114 120

8. 月刊 C MAGAZINE 2002年3月号

れほど問題はなく , 事実一般にも非常に広 と合わせて , 再び実数と虚数のペアに直し フェーズボコーダによるピッチスケ く用いられています。 ます。各正弦波でこの処理が終わったら , ーリング 最後に逆 F を行って , 周波数軸上のスペ 位相差計測法 こまでの話をまとめて , いよいよピッ クトル情報から時間軸上の情報に直して , さて , 「真の周波数」を求める方法はいろ チスケーリングの構成に移りましよう。今 信号として出力します。結局のところ , いろ考案されていますが , こでは位相差 回構成しようとするピッチスケーリングは , F ロ、で得られた各スペクトルの周波数情報 計測法と呼ばれる方法を用います。 FFT で のみをすり替えて再び逆 F で合成してい Fig. 33 のようなプロックダイアグラムとな 分解される各正弦波の真の周波数 ( 角速度 ) ります。まこと複雑な構成と思われるかも ることにほかなりません。 は Fig. 30 の複素平面上で原点を中心とした しれませんが , ソースはけっこうきれいに なお , スケーリング操作では注意すべき まとめられるので , ご安心ください。プロ べクトルの回転速度に相当します。位相差 ポイントがあります。自然界の信号はリニ 計測法ではこれをうまく利用してやるので ックダイアグラム上でデータの流れを順に アスケールではなく口グスケールで秩序立 す。 Fig. 32 を見てください。 追っていきましよう。 てられるので , スケーリング操作はあくま まず時間軸上のごく近い領域 ( お互いが まずは信号を F ドで 1 本 1 本の正弦波に分 でも掛け算で行ってください。もしも 1 オ 重なるくらい近く ) で 2 度 F ド I 、を行います。 解してしまいます。なお , F ドの結果は左 クタープアップなら , 1KHz を 2KHz に , 2K そして周波数を求めたい k 番目のスペクト 右対称なので N 本の正弦波すべてを処理す Hz を 4KHz にするということです。もしも ル同士の位相を比べるのです。この 2 つの るのは冗長です。よって , まず得られた N / これを足し算で IKHz を 2KHz に , 2KHz を 3 2 セットの ( R [ k ] , Ⅱ k ] ) という実数と虚数の KHz にするようなスケーリングをしてしま F ドフレームの時間差は図に示すように」 t ペアを式とアークタンジェントを使って うと , 元の波形では 2 倍音関係にあった IK で既知なので , この 2 つのフレーム間での k 振幅と位相に分解してしまいます。振幅情 番目のスペクトルの位相差はが」の時間 Hz と 2KHz のスペクトルが , シフト後は IK Hz と 3KHz という , 3 倍音関係に変更され でどれだけ位相が進むかによって理論値を 報はそのまま保存しておきます。位相情報 求めることができます。ところが実際に位 は次の位相差計測法で真の周波数に変換さ てしまいます。これは音楽信号で考えると , 相の差を計測してみると , たいていこの理 れます。そして次にシフトしたいぶんだけ , せつかくキレイな和音関係にあったものが 論値からずれるのです。実はこの「ずれ」 この真の周波数をシフトするスケーリング めちゃくちゃな不協和音になってしまうこ のぶんが真の周波数と / k とのずれにほかな 操作を行います。周波数がスケーリングで とを意味し , 調和関係の破壊につながりま らず , ここから真の周波数が推定できると きたら , 次にこの周波数を位相差計測法と す。 はまったく逆の処理を行い , 真の周波数が ちなみに , Fig. 23 のように F ド I 、などのフ いう仕組みです。 イルタバンクを使って各正弦波に分解し , なおこの」 t は短いほどよい結果が得ら FFT のグリッド上にうまく並ぶようにしま す。そして先ほど保存しておいた振幅情報 振幅・位相情報に分解し再合成していくよ れます。 正弦波 1 本 1 本の処理 振幅情報 実数部・虚数部に 分解 位相情報 Fig. 33 ピッチスケーリングの構成 正 の 本 位相情報 位相差計測法による 真の周波数測定 1 本の正弦波ー 周波数軸上での シフト [ = : 〕 ) = ・位相計測法の 周波数逆処理 周波数 情報 情報 同上 同上 逆 F 阿 - ー信号 各正弦波から合成 信号 各正弦波に分解 ソフトウェアによるオーティオデジタル信号処理のレシピ 31 特集 1 デジタルフィルタ

9. 月刊 C MAGAZINE 2002年3月号

作業です。 このあたりは , 単に時間領域で切り取っ たフレームの連結間隔で操作するという , ちょっと強引な OLA 系の方法の原理的な 問題といえるかもしれません。 アライメントボイントの問題 次にアライメントボイントについて詳し く見ていきます。 SOLA は先ほど解説した とおり , フレーム連結を間隔 S で行わず , 波形のつじつまの合う位置に多少ずらして 処理するものでした。しかしこの処理によ ってかえって OLA よりも悪化してしまう部 デジタルフィルタの 開発サイクル 筆者がデジタルフィルタを開発するとき の様子を少し紹介したいと思います。 デジタルフィルタの開発ツールには有名 な MATLAB があり , これを使うのが一般的 なようです。しかし筆者は MATLAB は使っ たことがありません。では , どのようにし て開発しているのでしようか ? ①設計するフィルタの伝達関数を資料か ら持ってくる ② Mathematica で理論計算を行いながら バラメータをいじってみて , バラメー タと周波数特性・位相特性の間の関係 をつかむ ③ VisuaI C + + などを用いて , 簡単なダイ アログペースアプリで実際に音を聞い てみる ④問題がなければ , お目当てのアプリケ ーションに組み込む といった開発サイクルで行っています。必 ずしも効率のよい方法とはいえませんが , 開発工程にブラックポックスの部分がなく , すべてひとつひとつ納得しながら実装まで いけるのがメリットだと思っています。 一般に Mathematica は遅いといわれてい ますが , デジタルフィルタの理論計算程度 ならもともとたいした演算量でもなく , と くに問題も起こりません。いかな周波数特 性といえども手計算や自分のプログラムで 数値計算をする気にはならないので , むし ろ Ma 物 ematica をたいへん重宝しています。 様子を簡単にプロットできるのもたいへん 便利な機能です。 Mathematica をお持ちの方は , 一度試し てみてください。 分もあります。 SO はの最大の欠点は , フ レーム間隔が必ずしも S にならないため , 式で定義される所望の変換比が厳密には 保証されないことです。 実はこのことはステレオファイルでは大 きな問題になってしまいます。なぜなら左 右チャンネル間で毎回のフレーム連結のア ライメントボイント lc 。が同じ値になるとは 一般には考えにくく , トータルで見た場合 にはチャンネル間で微妙にデータの長さが 違ってしまうことになります。実際のとこ ろ , ステレオファイルでは左右の位相情報 がけっこうシビアで , チャンネルごとに違 う km をとっていたのでは , ステレオ感も大 きく変化して ( というより破壊 ) しまいます。 このためステレオファイルの場合はどちら かのチャンネルでん , を求め , 他方もこのん , , でアラインするしかないようです。この際 lc, 。を評価しないほうのチャンネルではフレ ーム同士がベストボイントでつながらなく なり , 当然不利になります。しかし試聴を 重ねた結果 , それでもステレオ感がなくな るよりもずっとマシであるという結論に至 っています。なお , より低音が多く入って いるほうのチャンネルで k, 。 , を求めるのがよ EM-TSM の実装 いようです。 という内部キャッシュの末尾に連結されま nput は , 毎回確実に m—pLeftX , m—pRightX 入力バッフアである pLeftInput , pRightI を使っています リングクラス同様 , 内部でキャッシュ機構 CTSMUnit クラスも PART2 のピッチスケー とめてしか処理ができません。このため , フレームを使うので , データをある程度ま して設定する仕様になっています。爲 M は ネルの入力・出力バッフアをそれぞれ独立 この関数はステレオ専用です。左右チャン CT1meScaleUnit::Opearate 関数です。なお , この中で実際に時間軸を変更する関数は , ScaIeUnit クラスで行うようにしています。 アプリケーションでは , TSM をこの CTime 行います。 List 5 を見てください。サンプル ではいよいよ EM - TSM のソースの解説を コフ の超高速 信号処理バッケージ 本稿の各種のフィルタや計算ルーチンは , すべて「手」で書いています。しかし一般に は , 信号処理専用のライブラリ , フィルタ 演算部分などのライブラリを使うのが普通 です。その中で紹介したいのがコ n 回 SPL (Signal Processing Library) という信号処 理ライブラリ集です。この回 SPL は , 内 部で CPU を自動判定し , PentiumII, Ⅲ , 4 とどんどん拡張してきたマルチメディア 系の命令をフルに活用するようなコードに なっているようです。事実 , 各種のフィル タなどを演算させるしビックリするくら いの速度で動作します。 FFT の計算などに いたっては , 「本当に計算してるの ? 」と心 配になるようなパフォーマンスを見せてく れます。ライブラリはかなり実用的な内容 で物信号処理に必要な関数はたいていそろ っています。 特筆すべきは , DLL のカスタムビルドが 可能な点です。こうしたライブラリの場合 , 大きな DLL を数個アプリケーションに添え て配布しなくてはならないのが普通です。 しかしこの SPL ではライブラリから自分の 必要な関数のみをカスタム DLL として工ク スポートできるため , ライブラリ全体の大 きな DLL を配布する必要がありません。 れは再配布を考えたときにはメリットにな るでしよう。日本語のサイトもなぐ n 回リ リースのわりにはマイナーな存在のようで すが、ドキュメントも比較的しつかりして おり , きちんと使えば大きな武器になると 思います。ただ , 使いこなすにはある程度 の基本的な信号処理の知識は必要です。 URL は , http://www.intel.com/softwa 「 e/ products/perflib/spl/index. htm です。 す。そしてこの内部キャッシュのサイズが 1 フレーム以上処理できる大きさになると , いよいよ M を開始します。 TSM の最初で , まずべストアライメント ポイント k 。 , のサーチを式を用いて行いま す。このとき , km は Fig. 40 のように kMin と kMax ではさまれた領域内で見つけるよ うになっています。デフォルトでは士 400 ( マクロを使って , RANGE 値が定義されて いる ) としており , およそ 55Hz 程度の信号 の周期にあたります。この範囲は長ければ 長いほど低音の再現性がよくなりますが , 38 C MAGAZINE 2002 3

10. 月刊 C MAGAZINE 2002年3月号

コフ 6 電気工作からソフトでの信号処理へ 昔でいう「ラジオ少年」を , 真空管でなく トランジスタに置き換えてやっていたのが 筆者です。振り返ってみると , 小さいころ に父親 ( 真空管ラジオ少年だったらしい ) に見てもらいながらさんざん遊んだラジオ やオーディオの電気工作の知識が , 本稿の ような分野を扱ううえで少なからず役立っ ているように思えます。ハンダゴテを握り , 電気部品を組み合わせて昔はハードでやっ ていたことを , 今はそっくりそのままソフ トで「工作」していることになるのでしよう。 もっともテスターは Mathematica に変わ り , ハンダゴテは Visual C + + になったため , 火傷の心配はなくなりましたが。 そうした状況やツールは変わっても , 小 さいころ工作していたさまざまな経験がデ ジタル処理でも「実体験に基づく現象の把 握」というかたちで予備知識として役に立 っていることが少なくありません。 IIR フィ 1 月号拙稿 ) といった二種の神器 ( ? ) とでもいうような 各次元での操作方法すべてが手に入りまし た。あとは応用しだいです。想像力をはた らかせて , ぜひみなさんなりの用途を考え てみてください。こうしたマニアック ( ? ) な要素技術がどんなものに応用できるかを 考えているだけでも , とても楽しいもので す。 本特集では PC 上でのソフトウェアによ るオーディオアプリケーションを念頭に置 き , すぐに応用ができそうな各種の信号処 理技術を見てきました。信号処理というと 最近は MP3 などの圧縮系の技術に目がいき がちですが , 今回特集したような基礎技術 もたいへん重要なもので , 本格的なオーデ イオアプリケーションを作るためには欠か すことができないものです。もちろん今回 紹介したものはその中でごく一部にすぎず , オーディオ向けのデジタル信号処理にはま だまだ多くのものがあります。今回紹介し きれなかった部分についてはいずれ別の機 会に譲りたいと思います。 最後に ルタなどはその典型で , 電気部品でやって いたことをそっくりそのままソフトウェア に置き換えただけです。アナログ回路の恐 るべき現象 , 「発振」ですら起こるのですか ら。アナログ回路の体験なくいきなりデジ タル信号処理から入ると , こんな現象はわ けのわからない怪奇現象に映るかもしれま せん。しかしアナログ回路の経験があると 「発振しているな。過剰にフィードバック がかかっているのかもしれない。ちょっと 理論計算してみよう」という発想が生まれ , トラブルにすばやく対応できます。何しろ , アナログ回路で聞き慣れた「あの音」が PC のスピーカから聞こえるのですから。 こうした経験から , これからデジタル信 号処理をしようと考える人は , オーディオ 系の電気工作にもぜひチャレンジしていた だきたいと思うのです。きっと新しい発見 があるに違いありません。 以前はデジタル信号処理は専用 DSP ( Di gital Signal Processor) の独壇場だったので すが , 近年のマイクロプロセッサ技術のお かげで , 家庭用の PC でもソフトウェアで気 楽に実装できるようになってきました。添 付のサンプルアプリケーションでぜひ実際 に体験してみてください。 日本では PC 上のオーディオアプリケー ションという観点からの文献や資料は海外 に比べるとまだまだ少なく , 寂しい思いが します。今回の特集が読者のみなさんにこ うした分野への興味を持ってもらうきっか けとなれば幸いです。 ・参考文献 ・全般 Digital Signal Processing Second Editi0 n ,Sanjit K. Mitra , McGraw-HiII Higher Education ・「ディジタル信号処理の基礎』 , 三上直樹 著 , CQ 出版 ・ PARTI ・「ディジタルフィルタと信号処理』 , 谷萩 隆嗣著 , コロナ社 ・「 C 用ディジタル信号処理』 , 小畑秀文 / 幹 テフ材尠 康著 , コロナ社 ・「ソフトウェアで実現するディジタルグ ラフィックイコライザの設計と実装」 (TI nterface 』 2000 年 12 月号 ) , 小川要著 , CQ 出版 ℃ 00kb00k formulae for audio EQ biquad filter coefficients" (http://www.harmony- central.com/Computer/Programming/A udio-EQ-Cookbook. txt) , Robert Bristow- Johnson 著 A Digital Signal Processing Primer", Ke Ⅱ Steiglitz 著 , Addison-WesIey Publishing ・ PART2 ・『デジタルフーリエ変換ビギナーズ』 , 中村尚五著 , 東京電機大学出版局 ・「 FFT ( 高速フーリエ・コサイン・サイン 変換 ) の概略と設計法」 (http://momonga. t. u-tokyo. ac.jp/-ooura/fftman/index. html), 大浦拓哉著 'Stephan M. Sprenger's Audio DSP Page s" (http://www.dspdimension.com/) , S tephan M. Sprenger ・ PA RT3 ・「サンプリングレート変換によるスロー / 早送りアプリケーションの理論と実装」 ()C MAGAZINE 』 2002 年 1 月号 ) , 小川要 著 , ソフトバンクパプリッシング High Quality Tme Scale Modification fo r Speech" (Proc IEEE lnt. Conf. Acoustic, Speech, Signal processing, VOI. 1 (pp. 493 -496 ) , March 1985 ) , S. Roucos, AM. Wi lgus 著 Fast T1me ScaIe Modification using Enve lope-Matching Technique (EM-ISM) , " (P roc. of 1998 IEEE lnt. Sym. on Circuits & Systems(ISCAS) , vol. 5 (pp. 55 553 ) , M onterey, CA, USA, June 1998. ) , J. 、 v. C. W ・ ong, 0. C. Au, P. H. 、 v. Wong ・そのほか有用な www リソース Harmony Central@-Audi0 Programmin g , http://www.harmony-central.com/C omputer/Programming/ ・ music-dsp , 信号処理関連 ML ( 英語 ) , http://shoko.calarts.edu/%7Eglmrboy/m usicdsp/music-dsp. h加11 特集 1 デジタルフィルタソフトウェアによるオーティオデジタル信号処理のレシピ 41