第5章前半を追加

dsugawara/chapter05/q01.py

@@ -0,0 +1,34 @@
+import numpy as np
+
+
+def linear_array_vector(d, M, theta, f):
+    """直線状アレイのアレイマニフォールドベクトルを求める
+    Args:
+        d (double): アレイ間隔
+        M (int): マイク数
+        theta (int): 音源方向（y軸から反時計回りが正の向き）
+        f (int): 周波数
+    Return:
+        amv: アレイマニフォールドベクトル
+    """
+    c = 334.0
+    theta = np.pi / 2 + 2 * np.pi * (theta / 360)
+    u = np.array([np.sin(theta), np.cos(theta), 0]).T
+
+    p = np.zeros([M, 3])
+    for m in range(1, M + 1):
+        idx_m = m - 1
+        p[idx_m, 0] = ((m - 1) - (M - 1) / 2) * d
+    p = p.T
+    amv = np.empty([u.T.shape[0], p.shape[1]], dtype="complex")
+    amv = np.exp(2j * np.pi * f / c * np.dot(u.T, p))
+
+    return amv
+
+
+if __name__ == "__main__":
+    d = 0.05
+    M = 3
+    theta = 45
+    f = 1000
+    print(linear_array_vector(d, M, theta, f))

dsugawara/chapter05/q02.py

@@ -0,0 +1,35 @@
+import numpy as np
+
+
+def circular_array_vector(r, M, theta, f):
+    """円状アレイのアレイマニフォールドベクトルを求める
+    Args:
+        r (double): アレイ半径
+        M (int): マイク数
+        theta (int): 音源方向（y軸から反時計回りが正の向き）
+        f (int): 周波数
+    Return:
+        amv: アレイマニフォールドベクトル
+    """
+    c = 334.0
+    theta = np.pi / 2 + 2 * np.pi * (theta / 360)
+    u = np.array([np.sin(theta), np.cos(theta), 0]).T
+
+    p = np.zeros([M, 3])
+    for m in range(1, M + 1):
+        idx_m = m - 1
+        p[idx_m, 0] = r * np.sin(2 * np.pi / M * (m - 1))
+        p[idx_m, 1] = r * np.cos(2 * np.pi / M * (m - 1))
+    p = p.T
+    amv = np.empty([u.T.shape[0], p.shape[1]], dtype="complex")
+    amv = np.exp(2j * np.pi * f / c * np.dot(u.T, p))
+
+    return amv
+
+
+if __name__ == "__main__":
+    r = 0.05
+    M = 3
+    theta = 45
+    f = 1000
+    print(circular_array_vector(r, M, theta, f))

dsugawara/chapter05/q03.py

@@ -0,0 +1,54 @@
+import numpy as np
+
+
+def array_vector(crd, theta, f):
+    """アレイの座標からアレイマニフォールドベクトルを求める
+    Args:
+        crd (ndarray): アレイの座標
+        theta (int): 音源方向（y軸から反時計回りが正の向き）
+        f (int): 周波数
+    Return:
+        amv: アレイマニフォールドベクトル
+    """
+    M = crd.shape[0]
+    c = 334.0
+    theta = np.pi / 2 + 2 * np.pi * (theta / 360)
+    u = np.array([np.sin(theta), np.cos(theta), 0]).T
+
+    p = np.zeros([M, 3])
+    for m in range(1, M + 1):
+        idx_m = m - 1
+        r = crd[idx_m, 0] ** 2 + crd[idx_m, 1] ** 2
+        if r == 0:
+            sin, cos = 0, 0
+        else:
+            sin = crd[idx_m, 1] / r
+            cos = crd[idx_m, 0] / r
+        p[idx_m, 0] = r * sin
+        p[idx_m, 1] = r * cos
+    p = p.T
+    amv = np.empty([u.T.shape[0], p.shape[1]], dtype="complex")  # a[1, M]になるはず
+    amv = np.exp(2j * np.pi * f / c * np.dot(u.T, p))
+
+    return amv
+
+
+if __name__ == "__main__":
+    d = 0.05
+    theta = 45
+    f = 1000
+
+    # 直線状アレイの確認
+    liner = np.array([[-d, 0, 0], [0, 0, 0], [d, 0, 0]])
+    print(array_vector(liner, theta, f))
+
+    # 円状アレイの確認
+    r = 0.05
+    circular = np.array(
+        [
+            [r, 0, 0],
+            [r * np.cos(2 * np.pi / 3), r * np.sin(2 * np.pi / 3), 0],
+            [r * np.cos(2 * np.pi * 2 / 3), r * np.sin(2 * np.pi * 2 / 3), 0],
+        ]
+    )
+    print(array_vector(circular, theta, f))

dsugawara/chapter05/q04.py

@@ -0,0 +1,28 @@
+import numpy as np
+
+
+def scm(X):
+    """空間相関行列（spatial correation matrix）を求める
+    Args:
+        X (ndarray): 複素数行列
+    Return:
+        R: 空間相関行列
+    """
+    M, F, T = X.shape
+
+    R = np.empty([F, M, M], dtype="complex")
+    x = np.empty([F, M, T], dtype="complex")
+    for f in range(0, F):
+        for m in range(0, M):
+            x[f, m] = X[m, f]
+    for f in range(0, F):
+        R[f] = np.dot(x[f], np.conjugate(x[f].T)) / T
+    return R
+
+
+if __name__ == "__main__":
+    X1 = np.array([[1, -1j, -1, 1j], [2, -2j, -2, 2j], [3, -3j, -3, 3j]])
+    X2 = np.array([[4, -2j, 1, 0], [2, -1j, 0, 0], [1, -1j, 1, 0]])
+    X = np.array([X1, X2])
+
+    print(scm(X))

dsugawara/chapter05/q05.py

@@ -0,0 +1,103 @@
+import numpy as np
+
+
+def padding(L, S, x):
+    """N点の信号xにゼロ詰めを行う
+    Args:
+        L(int): 窓幅
+        S(int): シフト幅
+        x(ndarray): 入力信号
+    Return
+        x_pad: ゼロ詰め後の信号
+    """
+
+    x_pad = np.pad(x, (L - S, L - S))
+    temp = S - len(x_pad) % S
+    x_pad = np.pad(x_pad, (0, temp))
+    return x_pad
+
+
+def flame_div(L, S, x):
+    """N点の信号xをフレーム分割
+    Args:
+        L(int): 窓幅
+        S(int): シフト幅
+        x(ndarray): 入力信号
+    Return
+        x_div: フレーム分割したx
+    """
+    x_pad = padding(L, S, x)
+    N = len(x_pad)
+    T = ((N - L) // S) + 1
+    x_div = np.empty([T, L])
+    for t in range(0, T):
+        x_div[t] = x_pad[t * S : t * S + L]
+
+    return x_div
+
+
+def my_stft(L, S, x, win):
+    """短時間フーリエ変換
+    Args:
+        L(int): 窓幅
+        S(int): シフト幅
+        x(ndarray): 入力信号
+        w(ndarray): 窓関数
+    Return
+        X: 変換後の信号
+    """
+    x_div = flame_div(L, S, x)
+    T = x_div.shape[0]
+    X = np.empty([L // 2 + 1, T], dtype="complex")
+    for t in range(0, T):
+        X[:, t] = np.fft.rfft(x_div[t, :] * win)
+
+    return X
+
+
+def scm(X):
+    """空間相関行列（spatial correation matrix）を求める
+    Args:
+        X (ndarray): 複素数行列
+    Return:
+        R: 空間相関行列
+    """
+    M, F, T = X.shape
+
+    R = np.empty([F, M, M], dtype="complex")
+    x = np.empty([F, M, T], dtype="complex")
+    for f in range(0, F):
+        for m in range(0, M):
+            x[f, m] = X[m, f]
+    for f in range(0, F):
+        R[f] = np.dot(x[f], np.conjugate(x[f].T)) / T
+    return R
+
+
+if __name__ == "__main__":
+    X1 = np.array([[1, -1j, -1, 1j], [2, -2j, -2, 2j], [3, -3j, -3, 3j]])
+    X2 = np.array([[4, -2j, 1, 0], [2, -1j, 0, 0], [1, -1j, 1, 0]])
+    X = np.array([X1, X2])
+
+    print(scm(X))
+
+
+if __name__ == "__main__":
+
+    A = 1.0
+    fs = 16000
+    sec = 5
+
+    nA = A * np.random.randn(round(fs * sec))
+    nB = A * np.random.rand(round(fs * sec))
+
+    L = 512  # 窓長
+    S = 256  # シフト幅
+    win = np.hanning(L)  # 窓関数
+
+    nA_stft = my_stft(L, S, nA, win)
+    nB_stft = my_stft(L, S, nB, win)
+
+    R = scm(np.array([nA_stft, nB_stft]))
+
+    print(R[100].real)

dsugawara/chapter05/q06.py

@@ -0,0 +1,44 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def make_noise(signal, is_SNR):
+    """SNRを調整した雑音を生成する
+    Args:
+        signal(ndarray): 信号
+        is_SNR(double): SNR
+    Return:
+        ndarray: 雑音
+    """
+    noise = np.random.randn(len(signal))
+
+    noise = noise / np.sqrt(np.sum(noise**2))
+    noise = noise * np.sqrt(np.sum(signal**2))
+    noise = noise * 10 ** (-1 * is_SNR / 20)
+
+    return noise
+
+
+if __name__ == "__main__":
+    A = 1.0
+    fs = 16000
+    sec = 1
+    f = 440
+    t = np.arange(0, fs * sec) / fs
+
+    s = A * np.sin(2 * np.pi * f * t)
+
+    noise = make_noise(s, 10)
+
+    x1 = s + noise
+    s2 = np.pad(s, (10, 0))
+    x2 = s2[0 : len(s2) - 10] + noise
+    s3 = np.pad(s, (20, 0))
+    x3 = s3[0 : len(s3) - 20] + noise
+
+    plt.figure(figsize=[6.0, 4.0])
+    plt.plot(x1)
+    plt.plot(x2)
+    plt.plot(x3)
+    plt.xlim(0, fs * 0.01)
+    plt.savefig("q06.pdf")