1. No category

Scilab 演習
2015/5/13
以下の項目についてプログラミングし，ソース，プロットした図および結果
の考察をレポートにまとめなさい．図だけのレポートは評価しないので注意す
ること．
１． まず、以下の URL より音楽ソース sample.wav と、その音楽ソースに正
弦波の雑音が混入した wa2.wav をダウンロードしなさい。
http://dsp.cse.kyutech.ac.jp/lec/h27_dsp.html
wa2.wav ソースを聞くと正弦波が混入し聞きづらい音楽となっている。そこ
で、この正弦波の周波数をフーリエ変換で推定し、この周波数スペクトルを
除去するように FIR フィルタでフィルタリング（畳み込み）をしなさい。フ
ィルタリングする前後の波形と振幅スペクトルのプロットおよびスピーカ
（ヘッドフォン）で音を出して聴覚試験をしなさい。
レポートは、プログラム、実行結果（処理前後の時間波形、フィルタの
振幅特性、処理前後の信号の振幅スペクトル、音質の聴覚評価）および結果
の考察を加えなさい。
聴覚評価は、１（非常に悪い）－５（原音と同じ：非常に良い）までの５段
階で実施しなさい。
原音 sample.wav は聴覚試験の比較用だけであるので、プログラムには使用
しないこと。
２． 設計と解析の指針
－１ フィルタ設計
例えば 0.2Hz のカットオフ周波数を有する LPF は
[h]=wfir('lp',20,[0.2 0],'hm', [1 0]);
で設計できる。ここで’lp’はフィルタタイプであり、ここでは LPF（低域通
過フィルタ）を示す。 20 はインパルス応答長 N、[0.2 0]はカットオフ周波数
（２番目の要素 0 はデフォルト）であり、0.5 がナイキスト周波数（サンプリ
ング周波数の半分）を示す。’hm’および[1 0]は窓関数に関係する引数で今回は
デフォルト値とする。設計者は、カットオフ周波数を変えて、できるだけ正
弦波雑音を除去しつつ、良好な音質を維持するように設計すべきである。
可能ならば低域通過フィルタ LPF でなく、帯域阻止フィルタ BRF（BSF）
を使ってほしい。
－２ 畳み込みは、convol(x,h)で実行可能。
－３ フィルタの振幅特性の測定
インパルス応答ｈを fft コマンドでフーリエ変換し、絶対値(abs)を取る。
－４ 信号の振幅スペクトルの推定
フーリエ変換コマンド fft を利用するが、wa2.wav のデータサイズが多い
ため、取り込めない。そこで、例えば信号 x2 の最初の 100 点だけ取り込んで
fft する方法を以下に示しておく。
// 信号 x2 の振幅スペクトルを見る
for j = 1:1:100;
x22(j)=x2(j);
end
X2_1 = fft(x22,-1);
X2_2 = abs(X2_1);
// x 軸の設定
xlin = linspace(0,13463,100);
scf(1);
plot(xlin,X2_2);
// タイトル,x 軸ラベル,y 軸ラベル記述
xtitle("frequency", "frequency[Hz]", "amplitude");
なお、音楽ソースファイルの入出力方法は、別紙に記載しておく。
３． 設計条件
インパルス応答長 N は、以下のように学籍番号に応じた値とする。
N=(学籍番号の下一桁)＋10 とする．
例： 学籍番号 13230010 ⇒ 10
15230009 ⇒ 19
その他に設計者は、カットオフ周波数を適切に設定しなければならない。
注意点；
Scilab；
基本的に配布した資料で説明するコマンドで十分対応できる
演習であるが，必要な場合は N408 室に来室して質問するか，各自ヘル
プを参照のこと．
連絡先：inoue@dsp.cse.kyutech.ac.jp
提出日；5 月 27 日（水）授業中に提出。
レポートには，学科，学籍番号，氏名を書くこと．
別紙 Scilab での wav ファイルの扱い方
まず、Scilab のカレントディレクトリに処理する wav ファイルがあるかを確認する。
・ ロードの仕方
[x,y] = loadwave(‘ファイル 1.wav’);
または
x = loadwave(‘ファイル 1.wav’);
ファイル１は読み込む wav データ名
・ セーブの仕方
savewave(‘ファイル 2.wav’ ,x [, rate ]);
※ファイル２はファイル１とは別の名前にしないと、ファイル１.wav が上書
きされてしまうので注意。
※[, rate]は、wa2.wav のサンプリングレートが 13463Hz なので、その値を引数と
する。
// リスト例(このリストで雑音除去が出来るとは限らない！)：
[x,y] = loadwave('wa2.wav');
[h]=wfir('hp',20,[0.2 0],'hm', [0 1]);
[xm,fr] = frmag(h,1000);
scf(1);
plot(fr, xm)
[y1]= convol(h, x);
savewave('way.wav', y1, 13463);