Skip to content

Latest commit

 

History

History
151 lines (133 loc) · 24.7 KB

exercises.md

File metadata and controls

151 lines (133 loc) · 24.7 KB
Raw

注意: 「実装せよ」など特に明記されている場合を除いて，適当なライブラリを使用して構いません． 推奨ライブラリは下記の通りです．また，下記ライブラリの使用を前提とした問題も一部出題されます．

第 1 章: 基礎

  1. 正弦波の生成: 振幅 1, 周波数 440 Hz の正弦波をサンプリング周波数 16000 Hz で 3 秒分作成しプロットせよ．
  2. WAV ファイルの作成（モノラル）: 1.で作成した正弦波を 16bit PCM フォーマットで wav ファイルとして保存せよ．
  3. WAV ファイルの作成（ステレオ）: 振幅 1, 周波数 660 Hz の正弦波をサンプリング周波数 16000 Hz で 3 秒分作成し，1.で作成した信号と合わせて 2ch の wav ファイルとして保存せよ．
  4. 白色雑音の生成: ホワイトノイズをサンプリング周波数 16000 Hz で 3 秒分作成しプロットせよ．
  5. 信号の混合: 1.で作成した正弦波と 4.で作成したホワイトノイズと正弦波を混合してプロットせよ．
  6. SN 比: 信号長の等しい 2 個の信号 $s[n], x[n]; (n = 0, \dots, N-1), $ の信号対雑音比 (SN比) $10 \log \frac{\sum _{n=0} ^{N-1} s[n] ^2}{\sum _{n=0} ^{N-1} x[n] ^2}$ を計算する関数を実装せよ．
  7. SN 比を指定した信号の混合（実装）: SN 比を任意の値に設定できるようにホワイトノイズの振幅を調整する関数を実装せよ．元信号と所望の SN 比を入力として受け取り，ホワイトノイズを重畳した信号を出力すること．
  8. SN 比を指定した信号の混合（確認）: 8.で実装した関数を用いて，6.と同様にホワイトノイズと正弦波の混合信号を作成し wav ファイルとして保存せよ．ただし，SN 比が 6dB となるようにすること．
  9. 間引きによるダウンサンプリング: 8.で保存した wav ファイルを読み込み，サンプリング周波数を 8kHz に変換して保存せよ．
  10. 簡単なフィルタ処理: 9.の信号に対して 5 点移動平均フィルタを適用した結果と元の信号をプロットせよ．

第 2 章: スペクトル解析

  1. DFT/IDFT の実装: $N$ 点の信号 $x[n] , (n = 0, \dots , N-1)$ の離散フーリエ変換 (DFT) $$X[k] = \sum _{n=0} ^{N-1} x[n] \exp \left( -j \frac{2\pi k n}{N} \right) , (k = 0, \dots, N-1)$$ とその逆変換(IDFT)$$x[n] = \frac{1}{N} \sum _{k=0} ^{N-1} X[k] \exp \left( j \frac{2\pi k n}{N} \right) , (n = 0, \dots, N-1)$$を計算する関数を実装せよ．
  2. DFT の確認: 1.で実装した関数を用いて 8 点の単位インパルス信号 $\delta [n] = [1, 0, 0, 0, 0, 0, 0, 0]$ の DFT を計算せよ．
  3. IDFT の確認: 2.の結果の IDFT を計算しプロットせよ．
  4. スペクトルの確認: 2.の結果の振幅スペクトルおよび位相スペクトル $|X[k]|, \angle X[k] , (k = 0, \dots, N-1)$ をプロットせよ．
  5. 既存の実装との比較: 8 点の単位インパルス信号の DFT を numpy.fft.fft 関数を用いて計算し，2.の結果と比較せよ．なお，これ以降 DFT を計算する場合は numpy.fft.fft関数を使用してよい．
  6. 正弦波のスペクトル: 第 1 章 1.で作成した信号の DFT を計算し，振幅スペクトルと位相スペクトルをプロットせよ．（プロットする際はデシベル表記にするために 20 * log10(np.abs(X)) などのようにしてデシベル表記になおしてから計算すると見やすいです．）
  7. 窓関数: 次式で定義される $N$点の Hamming 窓を作成する関数を実装せよ． $$w[n] = 0.54 - 0.46 \cos (2\pi \frac{n}{N - 1}) , (n = 0, \dots, N-1)$$
  8. 窓関数の確認: 第 1 章 1.で作成した信号に対して 6.の窓関数を適用した信号をプロットせよ．
  9. 窓関数の特性: 第 1 章 1.で作成した信号と同じ長さの Hamming 窓を作成し，DFT を計算せよ．
  10. 窓関数とスペクトルの関係: 第 1 章 1.で作成した信号の DFT を $X[k] (k = 0, \dots, N-1)$，9.の結果を $Y[k] (k = 0, \dots, N-1)$ とする． $X[k]$ と $Y[k]$ の巡回畳み込み $Z[k] = \sum _{n=0} ^{N-1} X[n] Y[k-n] (k = 0, \dots, N-1)$ の IDFT を計算し 8.の結果と比較せよ．

第 3 章: 畳み込みと簡単なフィルタ処理

  1. 線形畳み込み: $N$ 点の信号 $x[n], h[n] (n = 0, \dots, N-1)$ の線形畳み込み $z[n] = \sum _{k = 0} ^{N-1} x[k] h[n - k], (n = 0, \dots, 2N - 2)$ を計算する関数を実装せよ．ただし， $n - k < 0$ または $n - k > N-1$ では信号の値は 0 とすること．
  2. 巡回畳み込み: $N$ 点の信号 $x[n], h[n] (n = 0, \dots, N-1)$ の巡回畳み込み $z[n] = \sum _{k = 0} ^{N-1} x[k] h[(n - k) \mod N], (n = 0, \dots, N-1)$ を計算する関数を実装せよ．ただし， $a \mod b$ は $a$ を $b$ で割った余りを表す．
  3. 零詰め＋巡回畳み込み: $N$ 点の信号 $x[n], y[n]$ に適切な零詰めを行って巡回畳み込みを計算することで，線形畳み込みと同じ結果を与える関数を実装せよ．
  4. 各種畳み込みの関係: 1., 2., 3. で実装した関数を用いて $x[n] = [4, 3, 2, 1]$ と $y[n] = [1, 0, -1, 0]$ の線形畳み込み，巡回畳み込み，零詰めを行った巡回畳み込みを計算し結果をプロットせよ．
  5. 差分方程式（再帰なし）: 次の差分方程式 $y[n] = 0.2 x[n] + 0.2 x[n-1] + 0.2 x[n-2] + 0.2 x[n-3] + 0.2 x[n-4]$ で表される信号 $y[n]$ をプロットせよ． ただし， $x[n]$ は単位インパルス信号として 10 点プロットせよ．
  6. 差分方程式（再帰あり）: 次の差分方程式 $y[n] = 0.3 y[n-1] + 0.4 x[n]$ で表される信号 $y[n]$ をプロットせよ． ただし， $x[n]$ は単位インパルス信号として 10 点プロットせよ．
  7. 差分方程式（一般系）: 一般の差分方程式 $a_0 y[n] = - \sum _{k = 1} ^{N} a_k y[n-k] + \sum _{k = 0} ^{M} b_k x[n - k]$ で表される信号を計算する関数を実装せよ．ただし，入力は $a_0, \dots, a_N (a_0 \neq 0)$ の配列， $b_0, \dots, b_M$ の配列，および入力信号 $x[n] (n = 1, \dots, L)$ とすること．
  8. 周波数応答: 次の差分方程式 $y[n] = - \sum _{k = 1} ^{N} a_k y[n-k] + \sum _{k = 0} ^{M} b_k x[n - k]$ で表される信号の周波数応答 $H(e^{j\omega}) = \frac{\sum _{k=0} ^M b_k e^{-j\omega k}}{1 + \sum _{k=1} ^{N} a_k e^{-j\omega k}}$ を計算する関数を実装せよ．ただし， $\omega$ はサンプリング周波数を $f_s$ としたときの正規化各周波数 $\omega = \frac{2 \pi f}{f_s} (0 \leq f \lt f_s)$ である．
  9. 周波数応答の確認（再帰なし）: 8.で実装した関数を用いて，5.の信号の周波数応答 $H(e^{j\omega})$ を計算し，振幅特性 $|H(e^{j\omega})|$ および位相特性 $\angle H(e^{j\omega})$ をプロットせよ．ただし，サンプリング周波数 $f_s = 16000$ とし，$f = 0, \frac{1}{N} f_s, \dots, \frac{N-1}{N} f_s$ まで計算せよ（$N$ は適当に決めてよい）．
  10. 周波数応答の確認（再帰あり）: 8.で実装した関数を用いて，6.の信号の周波数応答 $H(e^{j\omega})$ を計算し，振幅特性 $|H(e^{j\omega})|$ および位相特性 $\angle H(e^{j\omega})$ をプロットせよ．ただし，サンプリング周波数 $f_s = 16000$ とし，$f = 0, \frac{1}{N} f_s, \dots, \frac{N-1}{N} f_s$ まで計算せよ（$N$ は適当に決めてよい）．

第 4 章: 短時間フーリエ変換

  1. 零詰め: 窓幅 $L$ ，シフト幅 $S (\lt L)$，および $N$ 点の実数値信号 $x[n]$ $(n = 0, \dots, N-1)$ を入力とし，下記の手順で零詰めした信号を出力する関数を実装せよ．
    1. 入力の先頭および末尾に $L - S$ 個の零詰め: $[\underbrace{0, \dots, 0}{L - S}, x[0], \dots, x[N-1], \underbrace{0, \dots, 0}{L - S}]$
    2. 最終的な信号の長さが $S$ の倍数となるよう入力の末尾に零詰め
  2. フレーム分割: 窓幅 $L$ ，シフト幅 $S (\lt L)$，および $N$ 点の実数値信号 $x[n]$ $(n = 0, \dots, N-1)$ を入力とし，1.で実装した関数を用いて $L$ 点の信号（以下，フレームと呼ぶ）の列を出力する関数を実装せよ．
    1. $t$ $(t = 0, \dots, T-1)$ 番目のフレームを次で求める: $x_t [l] \coloneqq \tilde{x}[tS+l]$ $(l = 0, \dots, L-1)$, ただし， $\tilde{x}[n]$ は零詰め後の入力
  3. STFT の実装: 窓幅 $L$ ，シフト幅 $S (\lt L)$，窓関数 $w[l] (l = 0, \dots, L-1)$ ，および $N$ 点の実数値信号 $x[n]$ $(n = 0, \dots, N-1)$ を入力とし，下記の手順で $x[n]$ の短時間フーリエ変換を出力する関数を実装せよ．このとき，出力は $(\frac{L}{2}+1) \times T$ の複素数行列となることに注意せよ．
    1. 2.で実装した関数を用いて，各フレームに対して窓関数 $w[l] (l = 0, \dots, L-1)$ をかけ np.fft.rfft 関数を適用する．
  4. STFT の確認: 振幅 1, 周波数 440 Hz の正弦波をサンプリング周波数 16000 Hz で 0.1 秒分作成し，2.で実装した関数を用いてスペクトログラムを計算せよ．ただし，窓幅を 1000，窓関数を Hamming 窓，シフト幅を 500 とすること．さらに，振幅スペクトログラム・位相スペクトログラムをプロットせよ．プロットには matplotlib.pyplot.pcolormesh を使用すること．
  5. 合成窓: シフト幅 $S$ と $L$ 点の窓関数 $w[l] (l = 0, \dots, L-1)$ を入力とし，次の結果を出力する関数を実装せよ．なお，これは逆短時間フーリエ変換に用いる最適合成窓のひとつである．
    1. $Q \gets \frac{L}{S}$
    2. $w_\mathrm{s}[l] = \frac{w[l]}{\sum _{m = -(Q - 1)} ^{Q - 1} (w[l - m S]) ^ 2}$
  6. ISTFT の実装: シフト幅 $S$ および $F\times T$ の複素数行列 $X[f, t]$ $(f = 0, \dots, F-1,\quad t = 0, \dots, T-1)$ を入力とし，下記の手順で $M$ 点の信号 $\hat{x}[n]$ を出力する関数を実装せよ．
    1. 窓幅 $N$ を計算および出力信号の長さ $M$ を計算: $N \gets 2 (F - 1)$, $M \gets S (T - 1) + N$
    2. 出力信号を初期化: $\hat{x}[m] = 0 \quad \forall m$
    3. 各フレーム $t$ について逆 DFT を計算: $z_t[n] \gets \mathcal{F}^{-1}(X[0, t], \dots, X[F-1, t])[n]$ $(n = 0, \dots, N-1)$ ただし，右辺は $F$ 点の複素数値信号 $X[0, t], \dots, X[F-1, t]$ を逆 DFT した結果の $n$ 番目の要素を表す．逆 DFT には np.fft.irfft 関数を使用すること．これにより， $z[n]$ の点数は $N = 2(F-1)$ となることに注意せよ．
    4. Overlap add の計算: $\hat{x}[t S + n] = \hat{x}[t S + n] + w_\mathrm{s}[n] z_t[n]$ ただし， $w_\mathrm{s}[n]$ は 5.で計算される最適合成窓， $z_t[n]$は上で計算した逆 DFT である．
  7. ISTFT の確認: 6.で実装した関数を用いて，4.の結果の逆短時間フーリエ変換を計算し結果をプロットせよ．
  8. 合成窓の確認: 6.の実装において，合成窓をすべての $n$ に対して $w_\mathrm{s}[n] = 1$ とした場合の ISTFT を実装し，前問の結果と比較せよ．
  9. 不確定性原理の確認: 周波数 100 ~ 16000 Hz のチャープ信号をサンプリング周波数 16000 Hz で 1 秒分作成し，窓幅/シフト幅を 100/50, 200/100, 400/200, 800/400 と変えながらスペクトログラムをプロットし結果を比較せよ．なお，チャープ信号の作成には scipy.signal.chirp 関数を使用せよ．
  10. 時間周波数のインデクスと物理量との対応: スペクトログラム，サンプリング周波数，シフト幅を入力とし，3.の実装で得られるスペクトログラムの縦軸と横軸の単位をそれぞれ周波数と秒に変換する関数を実装せよ．これを用いて，4.の結果を再度プロットして確認せよ．

第 5 章: アレイ信号処理基礎

  1. アレイマニフォールドベクトル（直線状アレイ）: アレイ間隔 $d$ [m]，マイク数 $M$ ，音源方向 $\theta$，周波数 $f$ [Hz] を入力とし，2 次元平面における直線状等間隔アレイのアレイマニフォールドベクトル $\boldsymbol{a}$ を計算する関数を実装せよ．ただし， $c$ は音速を表し，本演習では $c = 334$ [m/s]とする．また， $\theta$ は $y$ 軸から反時計回りを正の向きにとるものとする．アレイマニフォールドベクトルの $m\quad (m=1,\dots,M)$ 番目の要素 $a_m$ は音源方向ベクトル $\boldsymbol{u}$ および$m$ 番目のマイクにおける位置ベクトル $\boldsymbol{p}m$ を用いて次式で計算される．これを用いて， $d = 0.05$ [m]，$M = 3, \theta = 45 ^{\circ}$, $f = 1000$ [Hz] の場合の結果を print 文で確認せよ． $$ \begin{align} a_m &= \exp \left( j \frac{2\pi f}{c} \boldsymbol{u} ^{\top} \boldsymbol{p}{m} \right) \ \boldsymbol{u} &= \begin{bmatrix} \sin \theta & \cos \theta & 0 \end{bmatrix} ^{\top} \ \boldsymbol{p}_m &= \begin{bmatrix} \left( (m - 1) - \frac{M-1}{2} \right) d & 0 & 0 \end{bmatrix} ^{\top} \end{align} $$
  2. アレイマニフォールドベクトル（円状アレイ）: 1.と同様に，アレイ半径 $r$ [m]，マイク数 $M$ ，音源方向 $\theta$，周波数 $f$ [Hz] を入力とし，2 次元平面における円状等間隔アレイのアレイマニフォールドベクトル $\boldsymbol{a}$ を計算する関数を実装せよ．これを用いて， $r = 0.05$ [m]，$M = 3, \theta = 45 ^{\circ}$, $f = 1000$ [Hz] の場合の結果を print 文で確認せよ． $$ \begin{align} \boldsymbol{p}_m &= \begin{bmatrix} r \sin (\frac{2\pi}{M}(m-1)) & r \cos (\frac{2\pi}{M}(m-1)) & 0 \end{bmatrix} ^{\top} \end{align} $$
  3. アレイマニフォールドベクトル（一般）: アレイの座標 $[x_1, y_1, 0], \dots, [x_M, y_M, 0]$ [m] ，音源方向 $\theta$，周波数 $f$ [Hz] を入力とし，2 次元平面における一般のマイクアレイのアレイマニフォールドベクトルを計算する関数を実装せよ．これを用いて，1.および 2.と同様の結果となるようなパラメータを与えて結果を print 文で確認せよ．
  4. 空間相関行列の計算: $M$ 個の $F\times T$ の複素数行列 $X_m[f, t]$ $(f = 0, \dots, F-1,\quad t = 0, \dots, T-1) \quad (m = 1,\dots, M)$ を入力とし，下記の手順で空間相関行列を計算する関数を実装せよ．出力は $F$ 個の $M\times M$ 複素数行列となることに注意すること．これを用いて， $X_1 = [[1, -1j, -1, 1j], [2, -2j, -2, 2j], [3, -3j, -3, 3j]]$, $X_2 = [[4, -2j, 1, 0], [2, -1j, 0, 0], [1, -j, 1, 0]]$ として結果を print 文で確認せよ． a. $\boldsymbol{x}{ft} \coloneqq \begin{bmatrix} X_1[f, t] & \dots & X_M[f, t] \end{bmatrix} ^{\top}$ b. $\boldsymbol{R}{f} \gets \frac{1}{T} \sum {t = 0} ^{T-1} \boldsymbol{x}{ft} \boldsymbol{x}_{ft} ^{\mathsf{H}} \quad (f = 0, \dots, F-1)$ （$\boldsymbol{x}^\mathsf{H}$ は複素数ベクトル $\boldsymbol{x}$ の Hermite 転置を表す）
  5. 空間相関行列の確認: サンプリング周波数 16000 [Hz]，長さ 5 [s] のホワイトノイズを 2 チャネル作成し，窓幅 512・シフト幅 256・Hann 窓で STFT した信号に対して，4.で実装した関数を用いて空間相関行列を計算せよ．結果を $\boldsymbol{R}{f} \quad (f = 0, \dots, 513)$ とするとき， $\boldsymbol{R}{100}$ の実部を print 分で確認せよ．
  6. 簡単な多チャネル観測シミュレーション: サンプリング周波数 16000 [Hz]，長さ 1 [s] で振幅 1・周波数 440 [Hz]の正弦波 $s[n]$ と， $s[n]$ に対して SN 比 10dB になるよう振幅を調整したホワイトノイズ $\varepsilon[n]$ を作成せよ．さらに，$x_1[n] = s[n] + \varepsilon[n], x_2[n] = s[n-10] + \varepsilon[n], x_3[n] = s[n-20] + \varepsilon[n]$ として 3 チャネルの信号を作成し結果をプロットせよ．プロットの範囲は [0, 0.01] [s] など適当に拡大せよ．
  7. 遅延和ビームフォーマ: 次の手順で 6.の信号を強調した結果をプロットせよ． a. 6.の信号をチャネルごとに窓幅 1024・シフト幅 512・Hann 窓で STFT した結果を $X_1[f,t], X_2[f, t], X_3[f, t]$ とする．ただし， $F$ は周波数ビンの総数， $T$ はフレーム数である． a. $\boldsymbol{x}_{ft} \gets \begin{bmatrix} X_1[f, t] & X_2[f, t] & X_3[f, t]\end{bmatrix}^{\top}$ b. 各周波数 $f = \frac{f_s / 2}{F-1} 0, \dots, \frac{f_s / 2}{F-1}$ ごとにビームフォーマ $\boldsymbol{w}_f = \frac{1}{3} \begin{bmatrix} \exp (-j 2\pi f \tau _1) & \exp (-j 2\pi f \tau 2) & \exp (-j 2\pi f \tau 3) & \end{bmatrix}$ を計算する．ただし， $\tau_1 = 0, \tau_2 = \frac{10}{f_s}, \tau_3 = \frac{20}{f_s}$ である． d. 各周波数ごとに $Y[f,t] = \boldsymbol{w}{f}^{\mathsf{H}} \boldsymbol{x}{ft}$ として強調信号を求める． d. $Y[f,t]$ を逆 STFT する．
  8. ビームパターン: 周波数領域におけるビームフォーマのフィルタ $\boldsymbol{w}f \quad (f = 0, \dots, F-1)$ ，マイクアレイの座標 $\boldsymbol{p}m \quad (m = 1, \dots, M)$，サンプリング周波数 $f_s$ を入力とし，次の手順でビームパターンを描画する関数を実装せよ．これを用いて，1.の直線状アレイにおける遅延和ビームフォーマのビームパターンをプロットせよ． a. 角度 $\theta = 0, \dots, 360$ [deg] および周波数 $f = \frac{f_s / 2}{F-1} 0, \dots, \frac{f_s / 2}{F-1} F$ について 3.でアレイマニフォールドベクトル $\boldsymbol{a}{f}(\theta)$ を計算する． b. 角度 $\theta = 0, \dots, 360$ [deg] について $\Psi (f, \theta) = \boldsymbol{w}{f}^{\mathsf{H}} \boldsymbol{a}_{f} (\theta)$ を計算する． c. 横軸を角度，縦軸を周波数として $20 \log _{10} \lvert \Psi (f, \theta) \rvert$ をカラーバーでプロットする．
  9. 空間サンプリング定理: 空間的エイリアシングを避けるためにはマイク間隔 $d$ を $d \leq \frac{c}{2f}$ となるように選ぶ必要がある．ここで， $c$ [m/s] は音速，$f$ [Hz] は周波数である．8.で実装した関数を用いて，マイク間隔を 2 [cm], 5 [cm], 10 [cm] と変化させながら 1.で求めた直線状アレイのビームパターンをプロットし，結果を確認せよ．
  10. 空間スペクトル: $M$ チャネルの $N$ 点実数値信号 $z_{m}[n] \quad (m = 1, \dots, M)$ を入力とし，次の手順で空間スペクトルを計算する関数を実装せよ．これを用いて，マイク数3，マイク間隔0.05 mの直線状ビームフォーマおよび 6.で用いた信号を入力とし，周波数ビンインデクス $f=20, 21, \dots, 30$ ごとに結果を描画せよ．6.の信号は窓幅 1024・シフト幅 512・Hann 窓で STFT すること． a. 入力信号の STFT を $Z_m[f,t] \quad (f = 0, \dots, F-1,; t = 0, \dots, T-1,; m = 1,\dots, M)$ とする． b. 5.で実装した関数を用いて $Z_m[f,t]$ の空間相関行列を計算し $\boldsymbol{R}{z}$ とする． c. 角度 $\theta = 0, \dots, 360$ [deg] についてビームフォーマ $\boldsymbol{w}(\theta)$ を計算する． d. 角度 $\theta = 0, \dots, 360$ [deg] について $P(\theta) = \boldsymbol{w} ^{\mathsf{H}} (\theta) \boldsymbol{R}{z} \boldsymbol{w} (\theta)$ を計算する． e. 横軸を角度，縦軸を $20\log _{10} \lvert P(\theta) \rvert$ としてプロットする．

参考文献

書籍

解説記事

Web サイト