Merge pull request #29 from woodnx/skotsugi/chapter05

Add chapter05

Author: taishi-n

skotsugi/chapter05/q01.py

@@ -0,0 +1,22 @@
+import numpy as np
+
+
+def linear_array_manifold_vector(d: float, M: int, theta: float, f: float):
+    c = 334
+    a = np.zeros(M, dtype=complex)
+    u = np.array([np.sin(theta), np.cos(theta), 0])
+
+    for m in range(M):
+        p = np.array([(m - (M - 1) / 2) * d, 0, 0])
+        a[m] = np.exp(1j * 2 * np.pi * f / c * (u @ p))
+
+    return a
+
+
+if __name__ == "__main__":
+    d = 0.05
+    M = 3
+    theta = np.pi / 4
+    f = 1000
+    print(linear_array_manifold_vector(d, M, theta, f))
+    # [(0.7868536952034816-0.617139580925277j), (1+0j), (0.7868536952034816+0.617139580925277j)]

skotsugi/chapter05/q02.py

@@ -0,0 +1,23 @@
+import numpy as np
+
+
+def circular_array_manifold_vector(r: float, M: int, theta: float, f: float):
+    c = 334
+    a = np.zeros(M, dtype=complex)
+    u = np.array([np.sin(theta), np.cos(theta), 0])
+
+    for m in range(M):
+        phi = 2 * np.pi * m / M
+        p = np.array([r * np.sin(phi), r * np.cos(phi), 0])
+        a[m] = np.exp(1j * 2 * np.pi * f / c * (u @ p))
+
+    return a
+
+
+if __name__ == "__main__":
+    r = 0.05
+    M = 3
+    theta = np.pi / 4
+    f = 1000
+    print(circular_array_manifold_vector(r, M, theta, f))
+    # [0.7868537 +0.61713958j 0.97051349+0.2410468j  0.6148926 -0.78861086j]

skotsugi/chapter05/q03.py

@@ -0,0 +1,36 @@
+import numpy as np
+
+
+def array_manifold_vector(X: np.ndarray, theta: float, f: float):
+    c = 334
+
+    M, _ = X.shape
+    a = np.zeros(M, dtype=complex)
+    u = np.array([np.sin(theta), np.cos(theta), 0])
+
+    for m in range(M):
+        p = X[m]
+        a[m] = np.exp(1j * 2 * np.pi * f / c * u @ p)
+
+    return a
+
+
+if __name__ == "__main__":
+    theta = np.pi / 4
+    r = 0.05
+    f = 1000
+
+    X = np.array([[-r, 0, 0], [0, 0, 0], [r, 0, 0]])
+    print(array_manifold_vector(X, theta, f))
+    # [0.7868537-0.61713958j 1.       +0.j         0.7868537+0.61713958j]
+
+    M = 3
+    Y = np.zeros((M, 3))
+    for m in range(M):
+        phi = 2 * np.pi * m / M
+        s = r * np.sin(phi)
+        t = r * np.cos(phi)
+        Y[m] = [s, t, 0]
+
+    print(array_manifold_vector(Y, theta, f))
+    # [0.7868537 +0.61713958j 0.97051349+0.2410468j  0.6148926 -0.78861086j]

skotsugi/chapter05/q04.py

@@ -0,0 +1,35 @@
+import numpy as np
+
+
+def spatial_correlation_matrix(X_m: np.ndarray) -> np.ndarray:
+    M, F, T = X_m.shape
+    R = np.zeros((F, M, M), dtype=complex)
+
+    for f in range(F):
+        sum_x = np.zeros((M, M), dtype=complex)
+        for t in range(T):
+            x_ft = np.array([X_m[:, f, t]]).T
+            sum_x += x_ft @ np.conj(x_ft.T)
+
+        R[f] = sum_x / T
+
+    return np.array(R)
+
+
+if __name__ == "__main__":
+    X = np.array(
+        [
+            [[1, -1j, -1, 1j], [2, -2j, -2, 2j], [3, -3j, -3, 3j]],
+            [[4, -2j, 1, 0], [2, -1j, 0, 0], [1, -1j, 1, 0]],
+        ]
+    )
+
+    print(spatial_correlation_matrix(X))
+    # [[[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]]]

skotsugi/chapter05/q05.py

@@ -0,0 +1,20 @@
+import numpy as np
+from scipy.signal import stft
+from q04 import spatial_correlation_matrix
+
+fs = 16000
+sec = 5
+L = 512
+S = 256
+
+np.random.seed(1)
+noise = np.random.normal(0, 1, size=(2, sec * fs))
+
+f, t, X = stft(noise, fs, nperseg=L, noverlap=L - S)
+
+R = spatial_correlation_matrix(X)
+
+print(R[100].real)
+
+# [[3.19811768e-03 7.08151469e-06]
+#  [7.08151469e-06 3.29345106e-03]]

skotsugi/chapter05/q06.png

@@ -0,0 +1,38 @@
+import numpy as np
+
+fs = 16000
+sec = 1
+A = 1
+f = 440
+snr = 10
+
+# make sin wave
+t = np.arange(fs * sec) / fs
+signal = A * np.sin(2 * np.pi * f * t)
+signal_length = len(signal)
+
+# make noise
+signal_power = np.mean(signal**2)
+snr_linear: float = 10.0 ** (snr / 10.0)
+noise_power = signal_power / snr_linear
+noise = np.sqrt(noise_power) * np.random.normal(0, 1, size=signal_length)
+
+# signalを2つ結合する，
+signal_twice = np.concatenate([signal, signal])
+
+x = np.array(
+    [
+        signal + noise,
+        signal_twice[10 : signal_length + 10] + noise,
+        signal_twice[20 : signal_length + 20] + noise,
+    ]
+)
+
+if __name__ == "__main__":
+    import matplotlib.pyplot as plt
+
+    for x_i in x:
+        plt.plot(t, x_i)
+
+    plt.xlim([0, 0.01])
+    plt.savefig("./skotsugi/chapter05/q06.png")

skotsugi/chapter05/q06.py

@@ -0,0 +1,31 @@
+import numpy as np
+from scipy.signal import stft, istft
+from q06 import x, fs
+
+
+L = 1024
+S = 512
+fbin, tbin, X = stft(x, fs, window="hann", nperseg=L, noverlap=L - S)
+
+N, F, T = X.shape
+
+Y = np.zeros((F, T), dtype=complex)
+tau = np.array([0, 10, 20]) / fs
+
+for i in range(F):
+    f = i * fs / 2 / (F - 1)
+    w = np.exp(-1j * 2 * np.pi * f * tau) / 3
+
+    for t in range(T):
+        Y[i, t] = np.conjugate(w) @ X[:, i, t]
+
+t, y = istft(Y, fs, window="hann", nperseg=L, noverlap=L - S)
+
+
+if __name__ == "__main__":
+    import matplotlib.pyplot as plt
+
+    plt.plot(t, y)
+
+    plt.xlim([0, 0.01])
+    plt.savefig("./skotsugi/chapter05/q07.png")

skotsugi/chapter05/q07.png

@@ -0,0 +1,43 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import array_manifold_vector
+from q07 import F
+
+
+def beam_pattern(w: np.ndarray, p: np.ndarray, fs):
+    F, M = w.shape
+
+    angle = np.arange(360)
+    theta = angle * np.pi / 180.0
+    freq = np.arange(F) * fs / 2 / (F - 1)
+
+    Phi = np.zeros((F, 360), dtype=complex)
+
+    for i, f in enumerate(freq):
+        for j, th in enumerate(theta):
+            a = array_manifold_vector(p, th, f)
+            w_f = w[i]
+
+            Phi[i, j] = w_f.T @ a
+
+    A = 20 * np.log10(np.abs(Phi))
+    return angle, freq, A
+
+
+if __name__ == "__main__":
+    d = 0.05
+    fs = 16000
+    M = 3
+    c = 334
+
+    # 直線状アレイ
+    p = np.array([[-d, 0, 0], [0, 0, 0], [d, 0, 0]])
+
+    # ステアリングベクトル
+    w = np.zeros((F, M), dtype=complex)
+    for f in range(F):
+        a_w = array_manifold_vector(p, 0, f)
+        w[f, :] = np.conj(a_w) / M
+
+    plt.pcolormesh(*beam_pattern(w, p, fs))
+    plt.savefig("./skotsugi/chapter05/q08.png")

skotsugi/chapter05/q07.py

@@ -0,0 +1,25 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import array_manifold_vector
+from q08 import beam_pattern, F
+
+fs = 16000
+M = 3
+c = 334
+ds = [0.02, 0.05, 0.1]
+fig, ax = plt.subplots(1, 3)
+
+for i, d in enumerate(ds):
+    # ステアリングベクトル
+    w = np.zeros((F, M), dtype=complex)
+    for f in range(F):
+        a_w = array_manifold_vector(p, 0, f)
+        w[f, :] = np.conj(a_w) / M
+
+    # 直線状アレイ
+    p = np.array([[-d, 0, 0], [0, 0, 0], [d, 0, 0]])
+
+    ax[i].pcolormesh(*beam_pattern(w, p, fs))
+
+plt.tight_layout()
+plt.savefig("./skotsugi/chapter05/q09.png")

skotsugi/chapter05/q08.png

@@ -0,0 +1,40 @@
+import numpy as np
+from scipy.signal import stft
+from q01 import linear_array_manifold_vector
+from q04 import spatial_correlation_matrix
+
+
+def spatial_spectrum(z: np.ndarray, f: float, fs: float):
+    M, N = z.shape
+    d = 0.05
+    L = 1024
+    S = 512
+
+    fbin, _, Z = stft(z, fs, window="hann", nperseg=L, noverlap=L - S)
+    R = spatial_correlation_matrix(Z)
+
+    P = np.zeros(360, dtype=complex)
+
+    for deg in range(360):
+        theta = deg * np.pi / 180.0
+        a = np.array([linear_array_manifold_vector(d, M, theta, fbin[f])])
+        w = a.T / M
+
+        P[deg] = (np.conj(w).T @ R[f] @ w).item()
+
+    return P
+
+
+if __name__ == "__main__":
+    import matplotlib.pyplot as plt
+    from q06 import x, fs
+
+    degs = np.arange(360)
+
+    for k in range(11):
+        f = k + 20
+
+        P = spatial_spectrum(x, f, fs)
+        plt.plot(degs, 20 * np.log10(np.abs(P)))
+
+    plt.savefig("./skotsugi/chapter05/q10.png")