Merge pull request #28 from 153hashimoto/khashimoto/chapter05

５章を追加　橋本

Author: taishi-n

khashimoto/chapter05/q01.py

@@ -0,0 +1,22 @@
+# アレイマニフォールドベクトル（直線状アレイ）
+
+import numpy as np
+
+
+def array_manifold_vector1(d, M, theta, f):
+    c = 334
+    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 = np.array([((m - 1) - (M - 1) / 2) * d, 0, 0]).T
+        a[m - 1] = np.exp(1j * 2 * np.pi * f / c * np.dot(u, p_m))
+    return a
+
+
+d = 0.05
+M = 3
+theta = np.pi / 4
+f = 1000
+print(array_manifold_vector1(d, M, theta, f))
+
+# [0.7868537-0.61713958j 1.       +0.j         0.7868537+0.61713958j]

khashimoto/chapter05/q02.py

@@ -0,0 +1,28 @@
+# アレイマニフォールドベクトル（円状アレイ）
+
+import numpy as np
+
+
+def array_manifold_vector2(r, M, theta, f):
+    c = 334
+    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 = np.array(
+            [
+                r * np.sin(2 * np.pi / M * (m - 1)),
+                r * np.cos(2 * np.pi / M * (m - 1)),
+                0,
+            ]
+        ).T
+        a[m - 1] = np.exp(1j * 2 * np.pi * f / c * np.dot(u, p_m))
+    return a
+
+
+r = 0.05
+M = 3
+theta = np.pi / 4
+f = 1000
+print(array_manifold_vector2(r, M, theta, f))
+
+# [0.7868537 +0.61713958j 0.97051349+0.2410468j  0.6148926 -0.78861086j]

khashimoto/chapter05/q03.py

@@ -0,0 +1,37 @@
+# アレイマニフォールドベクトル（一般）
+
+import numpy as np
+
+
+def array_manifold_vector3(p, theta, f):
+    c = 334
+    M = len(p)
+    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].T
+        a[m - 1] = np.exp(1j * 2 * np.pi * f / c * np.dot(u, p_m))
+    return a
+
+
+# 1.
+p = np.array([[-0.05, 0, 0], [0, 0, 0], [0.05, 0, 0]])
+theta = np.pi / 4
+f = 1000
+print(array_manifold_vector3(p, theta, f))
+
+# [0.7868537-0.61713958j 1.       +0.j         0.7868537+0.61713958j]
+
+# 2.
+p = np.array(
+    [
+        [0, 0.05, 0],
+        [0.05 * (np.sqrt(3) / 2), 0.05 * (-1 / 2), 0],
+        [0.05 * (-np.sqrt(3) / 2), 0.05 * (-1 / 2), 0],
+    ]
+)
+theta = np.pi / 4
+f = 1000
+print(array_manifold_vector3(p, theta, f))
+
+# [0.7868537 +0.61713958j 0.97051349+0.2410468j  0.6148926 -0.78861086j]

khashimoto/chapter05/q04.py

@@ -0,0 +1,38 @@
+# 空間相関行列の計算
+
+import numpy as np
+
+
+def spatial_correlation_matrix(*X):
+    M = len(X)
+    F = X[0].shape[0]
+    T = X[0].shape[1]
+    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
+
+
+X_1 = np.array([[1, -1j, -1, 1j], [2, -2j, -2, 2j], [3, -3j, -3, 3j]])
+X_2 = np.array([[4, -2j, 1, 0], [2, -1j, 0, 0], [1, -1j, 1, 0]])
+R = spatial_correlation_matrix(X_1, X_2)
+print(R[0])
+print(R[1])
+print(R[2])
+
+# # [[1.  +0.j 1.25+0.j]
+#  [1.25+0.j 5.25+0.j]]
+# [[4.  +0.j 1.5 +0.j]
+#  [1.5 +0.j 1.25+0.j]]
+# [[9.  +0.j 0.75+0.j]
+#  [0.75+0.j 0.75+0.j]]

khashimoto/chapter05/q05.py

@@ -0,0 +1,54 @@
+# 空間相関行列の確認
+
+import numpy as np
+from q04 import spatial_correlation_matrix
+
+
+def zeropad(L, S, x):
+    zero1 = np.zeros(L - S)
+    x_pad = np.concatenate([zero1, x])
+    zero2 = x_pad.size % S
+    if zero2 != 0:
+        x_pad = np.concatenate([x_pad, np.zeros(S - zero2)])  # 2.
+    x_pad = np.concatenate([x_pad, np.zeros(L - S)])
+    return x_pad
+
+
+def frame(L, S, x):
+    x_pad = zeropad(L, S, x)
+    T = (x_pad.size - L) // S + 1  # フレーム数
+    x_t = np.zeros((T, L))
+    for t in range(T):
+        for l in range(L):
+            x_t[t][l] = x_pad[t * S + l]
+    return x_t
+
+
+def stft(L, S, w, x):
+    x_t = frame(L, S, x)
+    T = len(x_t)
+    x_stft = np.zeros((T, L // 2 + 1), dtype=complex)
+    for t in range(T):
+        x_stft[t] = np.fft.rfft(x_t[t] * w)
+    return x_stft.T
+
+
+fs = 16000
+ts = 5
+t = np.arange(ts * fs) / fs
+x = np.zeros((ts * fs, 2))
+x[:, 0] = np.random.normal(size=ts * fs)
+x[:, 1] = np.random.normal(size=ts * fs)
+
+L = 512
+S = 256
+w = np.hanning(L)
+
+x0_stft = stft(L, S, w, x[:, 0])
+x1_stft = stft(L, S, w, x[:, 1])
+
+R_f = spatial_correlation_matrix(x0_stft, x1_stft)
+print(R_f[100].real)
+
+# [[191.38231688   5.05200591]
+#  [  5.05200591 201.98183266]]

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)