５章後半を追加

Author: 153hashimoto

khashimoto/chapter05/q06.py

@@ -0,0 +1,35 @@
+# 簡単な多チャネル観測シミュレーション
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+fs = 16000
+sec = 1
+A = 1
+f = 440
+
+t = np.arange(fs * sec) / fs
+s = A * np.sin(2 * np.pi * f * t)  # 正弦波
+
+snr = 10
+noise = np.random.randn(fs * sec)
+noise /= np.sqrt(np.sum(noise**2))
+noise *= np.sqrt(np.sum(s**2))
+a = 10 ** (-snr / 20)
+noise = a * noise  # ホワイトノイズ
+
+x1 = np.zeros(fs * sec)
+x2 = np.zeros(fs * sec)
+x3 = np.zeros(fs * sec)
+for n in range(fs * sec):
+    x1[n] = s[n] + noise[n]
+    x2[n] = s[n - 10] + noise[n]
+    x3[n] = s[n - 20] + noise[n]
+
+
+plt.plot(t, x1)
+plt.plot(t, x2)
+plt.plot(t, x3)
+plt.xlim(0, 0.01)
+plt.show()

khashimoto/chapter05/q07.py

@@ -0,0 +1,73 @@
+# 遅延和ビームフォーマ
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import signal
+
+
+# 6.
+fs = 16000
+sec = 1
+A = 1
+f = 440
+
+ts = np.arange(fs * sec) / fs
+s = A * np.sin(2 * np.pi * f * ts)  # 正弦波
+
+snr = 10
+noise = np.random.randn(fs * sec)
+noise /= np.sqrt(np.sum(noise**2))
+noise *= np.sqrt(np.sum(s**2))
+a = 10 ** (-snr / 20)
+noise = a * noise  # ホワイトノイズ
+
+x1 = np.zeros(fs * sec)
+x2 = np.zeros(fs * sec)
+x3 = np.zeros(fs * sec)
+for n in range(fs * sec):
+    x1[n] = s[n] + noise[n]
+    x2[n] = s[n - 10] + noise[n]
+    x3[n] = s[n - 20] + noise[n]
+
+
+# 7.
+L = 1024
+S = 512
+window = np.hanning(L)
+# a.
+f_, t_, X1 = signal.stft(x1, fs, window, L, L - S)
+f_, t_, X2 = signal.stft(x2, fs, window, L, L - S)
+f_, t_, X3 = signal.stft(x3, fs, window, L, L - S)
+
+
+F, T = X1.shape
+f = np.arange(F) * (fs / 2) / (F - 1)
+tau1 = 0
+tau2 = 10 / fs
+tau3 = 20 / fs
+
+Y = np.zeros((F, T), dtype=complex)
+for i in range(F):
+    # c.
+    w_f = (
+        np.array(
+            [
+                np.exp(-1j * 2 * np.pi * f[i] * tau1),
+                np.exp(-1j * 2 * np.pi * f[i] * tau2),
+                np.exp(-1j * 2 * np.pi * f[i] * tau3),
+            ]
+        )
+        / 3
+    )
+    for j in range(T):
+        # b.
+        x_ft = np.array([X1[i, j], X2[i, j], X3[i, j]]).T
+        # d.
+        Y[i, j] = np.conj(w_f.T) @ x_ft
+
+# e.
+t_Y, Y_istft = signal.istft(Y, fs, window, L, L - S)
+
+plt.plot(t_Y, Y_istft)
+plt.xlim(0, 0.01)
+plt.show()

khashimoto/chapter05/q08.py

@@ -0,0 +1,64 @@
+# ビームパターン
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+# 3.
+def array_manifold_vector3(p, theta, f):
+    c = 334
+    M = p.shape[0]
+    a = np.zeros(M, dtype=complex)
+    u = np.array([np.sin(theta), np.cos(theta), 0]).T
+    for m in range(1, M + 1):
+        p_m = p[m - 1]
+        a[m - 1] = np.exp(1j * 2 * np.pi * f / c * (u.T @ p_m))
+    return a
+
+
+def beampattern(w_f, p_m, fs):
+    F = w_f.shape[0]
+    f = np.arange(F) * (fs / 2) / (F - 1)
+    psi = np.zeros((F, 361), dtype=complex)
+
+    for i in range(F):
+        for theta in range(361):
+            theta_rad = theta * np.pi / 180
+            # a.
+            a_f = array_manifold_vector3(p_m, theta_rad, f[i])
+            # b.
+            psi[i, theta] = np.conj(w_f[i, :].T) @ a_f
+
+    Theta = np.arange(361)
+    # c.
+    plt.pcolormesh(Theta, f, 20 * np.log10(np.abs(psi)))
+    plt.xlabel("Angle[°]")
+    plt.ylabel("Frequency[Hz]")
+    plt.show()
+
+
+c = 334
+fs = 16000
+L = 1024
+F = L // 2 + 1
+f = np.arange(F) * (fs / 2) / (F - 1)
+
+# 1. の直線状アレイにおける遅延和ビームフォーマのビームパターンをプロット
+d = 0.05
+M = 3
+theta = np.pi / 4
+tau1 = 0
+tau2 = d * np.sin(theta) / c
+tau3 = 2 * d * np.sin(theta) / c
+
+w_f = np.zeros((F, M), dtype=complex)
+for i in range(F):
+    w_f[i] = np.array(
+        [
+            np.exp(-1j * 2 * np.pi * f[i] * tau1),
+            np.exp(-1j * 2 * np.pi * f[i] * tau2),
+            np.exp(-1j * 2 * np.pi * f[i] * tau3),
+        ]
+    )
+p_m = np.array([[-d, 0, 0], [0, 0, 0], [d, 0, 0]])
+beampattern(w_f, p_m, fs)

khashimoto/chapter05/q09.py

@@ -0,0 +1,116 @@
+# 空間サンプリング定理
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+# 3.
+def array_manifold_vector3(p, theta, f):
+    c = 334
+    M = p.shape[0]
+    a = np.zeros(M, dtype=complex)
+    u = np.array([np.sin(theta), np.cos(theta), 0]).T
+    for m in range(1, M + 1):
+        p_m = p[m - 1]
+        a[m - 1] = np.exp(1j * 2 * np.pi * f / c * (u.T @ p_m))
+    return a
+
+
+# 8.
+def beampattern(w_f, p_m, fs):
+    F = w_f.shape[0]
+    f = np.arange(F) * (fs / 2) / (F - 1)
+    psi = np.zeros((F, 361), dtype=complex)
+
+    for i in range(F):
+        for theta in range(361):
+            theta_rad = theta * np.pi / 180
+            # a.
+            a_f = array_manifold_vector3(p_m, theta_rad, f[i])
+            # b.
+            psi[i, theta] = np.conj(w_f[i, :].T) @ a_f
+
+    return psi
+
+
+c = 334
+fs = 16000
+L = 1024
+F = L // 2 + 1
+f = np.arange(F) * (fs / 2) / (F - 1)
+M = 3
+theta = np.pi / 4
+
+
+# マイク間隔2[cm]
+d = 0.02
+tau1 = 0
+tau2 = d * np.sin(theta) / c
+tau3 = 2 * d * np.sin(theta) / c
+
+w_f = np.zeros((F, M), dtype=complex)
+for i in range(F):
+    w_f[i] = np.array(
+        [
+            np.exp(-1j * 2 * np.pi * f[i] * tau1),
+            np.exp(-1j * 2 * np.pi * f[i] * tau2),
+            np.exp(-1j * 2 * np.pi * f[i] * tau3),
+        ]
+    )
+p_m = np.array([[-d, 0, 0], [0, 0, 0], [d, 0, 0]])
+psi1 = beampattern(w_f, p_m, fs)
+
+# マイク間隔5[cm]
+d = 0.05
+tau1 = 0
+tau2 = d * np.sin(theta) / c
+tau3 = 2 * d * np.sin(theta) / c
+
+w_f = np.zeros((F, M), dtype=complex)
+for i in range(F):
+    w_f[i] = np.array(
+        [
+            np.exp(-1j * 2 * np.pi * f[i] * tau1),
+            np.exp(-1j * 2 * np.pi * f[i] * tau2),
+            np.exp(-1j * 2 * np.pi * f[i] * tau3),
+        ]
+    )
+p_m = np.array([[-d, 0, 0], [0, 0, 0], [d, 0, 0]])
+psi2 = beampattern(w_f, p_m, fs)
+
+# マイク間隔10[cm]
+d = 0.1
+tau1 = 0
+tau2 = d * np.sin(theta) / c
+tau3 = 2 * d * np.sin(theta) / c
+
+w_f = np.zeros((F, M), dtype=complex)
+for i in range(F):
+    w_f[i] = np.array(
+        [
+            np.exp(-1j * 2 * np.pi * f[i] * tau1),
+            np.exp(-1j * 2 * np.pi * f[i] * tau2),
+            np.exp(-1j * 2 * np.pi * f[i] * tau3),
+        ]
+    )
+p_m = np.array([[-d, 0, 0], [0, 0, 0], [d, 0, 0]])
+psi3 = beampattern(w_f, p_m, fs)
+
+Theta = np.arange(361)
+
+plt.subplot(1, 3, 1)
+plt.pcolormesh(Theta, f, 20 * np.log10(np.abs(psi1)))
+plt.xlabel("Angle[°]")
+plt.ylabel("Frequency[Hz]")
+
+plt.subplot(1, 3, 2)
+plt.pcolormesh(Theta, f, 20 * np.log10(np.abs(psi2)))
+plt.xlabel("Angle[°]")
+plt.ylabel("Frequency[Hz]")
+
+plt.subplot(1, 3, 3)
+plt.pcolormesh(Theta, f, 20 * np.log10(np.abs(psi3)))
+plt.xlabel("Angle[°]")
+plt.ylabel("Frequency[Hz]")
+
+plt.show()

khashimoto/chapter05/q10.py

@@ -0,0 +1,101 @@
+# 空間スペクトル
+
+import numpy as np
+from scipy import signal
+import matplotlib.pyplot as plt
+
+
+# 4.
+def spatial_correlation_matrix(X):
+    M, F, T = X.shape
+    R = np.zeros((F, M, M), dtype=complex)
+    for f in range(F):
+        sum = np.zeros((M, M), dtype=complex)
+        for t in range(T):
+            # a.
+            x_ft = np.zeros(M, dtype=complex)
+            for m in range(M):
+                x_ft[m] = X[m, f, t]
+            x_ft = np.array([x_ft]).T
+
+            # b.
+            sum = np.add(sum, x_ft @ np.conj(x_ft.T))
+        R[f] = sum / T
+    return R
+
+
+def spatialspec(z_m, M, d, f):
+    c = 334
+    fs = 16000
+    L = 1024
+    S = 512
+    window = np.hanning(L)
+
+    F = L // 2 + 1
+    T = (fs + L - S) // S + 1
+    # a.
+    Z_m = np.zeros((M, F, T), dtype=complex)
+    for m in range(M):
+        f_, t_, Z_m[m] = signal.stft(z_m[m], fs, window, L, L - S)
+
+    # b.
+    R_z = spatial_correlation_matrix(Z_m)
+
+    # c. d.
+    P = np.zeros(361, dtype=complex)
+    for theta in range(361):
+        theta_rad = theta * np.pi / 180
+        tau1 = 0
+        tau2 = d * np.sin(theta_rad) / c
+        tau3 = 2 * d * np.sin(theta_rad) / c
+        w = np.array(
+            [
+                np.exp(-1j * 2 * np.pi * f_[f] * tau1),
+                np.exp(-1j * 2 * np.pi * f_[f] * tau2),
+                np.exp(-1j * 2 * np.pi * f_[f] * tau3),
+            ]
+        )
+
+        P[theta] = np.conj(w.T) @ R_z[f] @ w
+    return P
+
+
+# 6.
+fs = 16000
+sec = 1
+A = 1
+f = 440
+
+ts = np.arange(fs * sec) / fs
+s = A * np.sin(2 * np.pi * f * ts)  # 正弦波
+
+snr = 10
+np.random.seed(0)
+noise = np.random.randn(fs * sec)
+noise /= np.sqrt(np.sum(noise**2))
+noise *= np.sqrt(np.sum(s**2))
+a = 10 ** (-snr / 20)
+noise = a * noise  # ホワイトノイズ
+
+x1 = np.zeros(fs * sec)
+x2 = np.zeros(fs * sec)
+x3 = np.zeros(fs * sec)
+for n in range(fs * sec):
+    x1[n] = s[n] + noise[n]
+    x2[n] = s[n - 10] + noise[n]
+    x3[n] = s[n - 20] + noise[n]
+
+
+# 10.
+z_m = np.array([x1, x2, x3])
+M = 3
+d = 0.05
+Theta = np.arange(361)
+
+freq = np.arange(20, 31)
+for i in range(11):
+    P = spatialspec(z_m, M, d, freq[i])
+    plt.subplot(6, 2, i + 1)
+    plt.plot(20 * np.log10(np.abs(P)))
+
+plt.show()