Merge pull request #26 from kuu8902/kkoiso/chapter05

５章前半を追加　小磯

Author: taishi-n

kkoiso/chapter05/q01.py

@@ -0,0 +1,21 @@
+import numpy as np
+
+
+def array_manifold_vector_linear(d, M, theta, f, c=334):
+    theta_rad = np.deg2rad(theta)
+    u = np.array([np.sin(theta_rad), np.cos(theta_rad), 0])
+    a = np.zeros(M, dtype=complex)
+    for m in range(M):
+        p_m = np.array([((m - 1) - (M - 1) / 2) * d, 0, 0])
+        a[m] = np.exp(1j * 2 * np.pi * f / c * np.dot(u, p_m))
+    return a
+
+
+d = 0.05
+M = 3
+theta = 45
+f = 1000
+
+
+a_linear = array_manifold_vector_linear(d, M, theta, f)
+print(a_linear)

kkoiso/chapter05/q02.py

@@ -0,0 +1,23 @@
+import numpy as np
+
+
+def array_manifold_vector_circular(r, M, theta, f, c=334):
+    theta_rad = np.deg2rad(theta)
+    u = np.array([np.sin(theta_rad), np.cos(theta_rad), 0])
+    a = np.zeros(M, dtype=complex)
+    for m in range(M):
+        p_m = np.array(
+            [r * np.sin(2 * np.pi * m / M), r * np.cos(2 * np.pi * m / M), 0]
+        )
+        a[m] = np.exp(1j * 2 * np.pi * f / c * np.dot(u, p_m))
+    return a
+
+
+r = 0.05
+M = 3
+theta = 45
+f = 1000
+
+
+a_circular = array_manifold_vector_circular(r, M, theta, f)
+print(a_circular)

kkoiso/chapter05/q03.py

@@ -0,0 +1,31 @@
+import numpy as np
+
+
+def array_manifold_vector_general(array_add, theta, f, c=334):
+    theta_rad = np.deg2rad(theta)
+    u = np.array([np.sin(theta_rad), np.cos(theta_rad), 0])
+    M = len(array_add)
+    a = np.zeros(M, dtype=complex)
+    for m in range(M):
+        p_m = np.array([array_add[m][0], array_add[m][1], 0])
+        a[m] = np.exp(1j * 2 * np.pi * f / c * np.dot(u, p_m))
+    return a
+
+
+d = 0.05
+M = 3
+theta = 45
+f = 1000
+linear_add = [((m - 1) * d, 0) for m in range(M)]
+
+r = 0.05
+circular_add = [
+    (r * np.sin(2 * np.pi * m / M), r * np.cos(2 * np.pi * m / M)) for m in range(M)
+]
+
+
+a_linear_general = array_manifold_vector_general(linear_add, theta, f)
+print(a_linear_general)
+
+a_circular_general = array_manifold_vector_general(circular_add, theta, f)
+print(a_circular_general)

kkoiso/chapter05/q04.py

@@ -0,0 +1,21 @@
+import numpy as np
+
+
+def spatial_correlation_matrix(X):
+    M, F, T = X.shape
+    R = np.zeros((F, M, M), dtype=complex)
+    for f in range(F):
+        for t in range(T):
+            x_ft = X[:, f, t]
+            R[f] += np.power(t, T - 1) * x_ft * np.conj(x_ft).T
+        R[f] /= T
+    return R
+
+
+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.stack((X1, X2), axis=0)
+
+
+R = spatial_correlation_matrix(X)
+print(R)

kkoiso/chapter05/q05.py

@@ -0,0 +1,20 @@
+import numpy as np
+from scipy.signal import stft
+import matplotlib.pyplot as plt
+from q04 import spatial_correlation_matrix
+
+fs = 16000
+T = 5
+N = fs * T
+np.random.seed(0)
+white_noise = np.random.normal(0, 1, (2, N))
+
+f, t, Zxx1 = stft(white_noise[0], fs, window="hann", nperseg=512, noverlap=256)
+f, t, Zxx2 = stft(white_noise[1], fs, window="hann", nperseg=512, noverlap=256)
+X = np.stack((Zxx1, Zxx2), axis=0)
+
+
+R = spatial_correlation_matrix(X)
+
+
+print(R[100].real)

kkoiso/chapter05/q06.py

@@ -0,0 +1,33 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+fs = 16000
+T = 1
+N = fs * T
+t = np.linspace(0, T, N, endpoint=False)
+f = 440
+a = 1
+
+s = a * np.sin(2 * np.pi * f * t)
+
+
+np.random.seed(0)
+noise = np.random.normal(0, 1, N)
+SNR = 10
+noise_amplitude = a / (10 ** (SNR / 20))
+epsilon = noise_amplitude * noise
+
+
+x1 = s + epsilon
+x2 = np.roll(s, 10) + epsilon
+x3 = np.roll(s, 20) + epsilon
+
+x = x1 + x2 + x3
+plt.figure(figsize=(10, 6))
+plt.plot(t[: int(0.01 * fs)], x[: int(0.01 * fs)])
+plt.xlabel("Time [s]")
+plt.ylabel("Amplitude")
+plt.title("Multichannel Signals")
+plt.legend()
+plt.show()

kkoiso/chapter05/q07.py

@@ -0,0 +1,54 @@
+from scipy.signal import istft, stft
+import numpy as np
+import matplotlib.pyplot as plt
+
+fs = 16000
+T = 1
+N = fs * T
+t = np.linspace(0, T, N, endpoint=False)
+f = 440
+a = 1
+
+s = a * np.sin(2 * np.pi * f * t)
+
+
+np.random.seed(0)
+noise = np.random.normal(0, 1, N)
+SNR = 10
+noise_amplitude = a / (10 ** (SNR / 20))
+epsilon = noise_amplitude * noise
+x = (s + np.roll(s, 10) + np.roll(s, 20)) / 3
+
+x1 = s + epsilon
+x2 = np.roll(s, 10) + epsilon
+x3 = np.roll(s, 20) + epsilon
+
+
+f, t, X1 = stft(x1, fs, window="hann", nperseg=1024, noverlap=512)
+f, t, X2 = stft(x2, fs, window="hann", nperseg=1024, noverlap=512)
+f, t, X3 = stft(x3, fs, window="hann", nperseg=1024, noverlap=512)
+
+X = np.stack((X1, X2, X3), axis=0)
+F, T = X1.shape
+
+
+tau = np.array([0, 10 / fs, 20 / fs])
+Y = np.zeros((F, T), dtype=complex)
+for f in range(F):
+    w_f = 1 / 3 * np.exp(-1j * 2 * np.pi * f * tau[:, None] * fs / 2 / F)
+    for t in range(T):
+        x_ft = X[:, f, t]
+        Y[f, t] = np.conj(w_f.T).dot(x_ft)
+
+
+_, y = istft(Y, fs, window="hann", nperseg=1024, noverlap=512)
+t = np.linspace(0, 1, 1 * fs, endpoint=False)
+
+
+plt.figure(figsize=(10, 6))
+plt.plot(t[: int(0.01 * fs)], y[: int(0.01 * fs)])
+plt.xlim(0, 0.01)
+plt.xlabel("Time [s]")
+plt.ylabel("Amplitude")
+plt.title("Enhanced Signal by Delay Sum Beamformer")
+plt.show()

kkoiso/chapter05/q08.py

@@ -0,0 +1,40 @@
+from scipy.signal import istft, stft
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import array_manifold_vector_general
+
+
+def plot_beam_pattern(w_f, p, fs):
+    fs_2 = fs / 2
+    angles = np.linspace(0, 360, 361)
+    F = w_f.shape[0]
+    Psi = np.zeros((F, len(angles)), dtype=complex)
+    for f in range(F):
+        for i, theta in enumerate(angles):
+            a_f_theta = array_manifold_vector_general(p, theta, f * (fs_2 / (F - 1)))
+            Psi[f, i] = np.conj(w_f[f]).dot(a_f_theta)
+    plt.figure(figsize=(10, 6))
+    plt.imshow(20 * np.log10(np.abs(Psi)), aspect="auto", extent=[0, 360, 0, fs / 2])
+    plt.colorbar(label="Amplitude [dB]")
+    plt.xlabel("Angle [degrees]")
+    plt.ylabel("Frequency [Hz]")
+    plt.title("Beam Pattern")
+    plt.show()
+
+
+d = 0.05
+M = 3
+theta = 45
+F = 1000
+linear_coords = [((m - 1) * d, 0) for m in range(M)]
+fs = 16000
+
+
+w_f = np.zeros((F, M), dtype=complex)
+for f in range(F):
+    fs_2 = fs / 2
+    tau = np.array([0, 10 / fs, 20 / fs])
+    w_f[f] = 1 / 3 * np.exp(-1j * 2 * np.pi * f * tau * fs_2 / (F - 1))
+
+
+plot_beam_pattern(w_f, linear_coords, fs)

kkoiso/chapter05/q09.py

@@ -0,0 +1,19 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q08 import plot_beam_pattern
+
+
+M = 3
+theta = 45
+F = 1000
+fs = 16000
+
+
+for d in [0.02, 0.05, 0.10]:
+    linear_coords = [((m - 1) * d, 0) for m in range(M)]
+    w_f = np.zeros((F, M), dtype=complex)
+    for f in range(F):
+        tau = np.array([0, 10 / fs, 20 / fs])
+        fs_2 = fs / 2
+        w_f[f] = 1 / 3 * np.exp(-1j * 2 * np.pi * f * tau * fs_2 / (F - 1))
+    plot_beam_pattern(w_f, linear_coords, fs)

kkoiso/chapter05/q10.py

@@ -0,0 +1,55 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.signal import stft
+from q04 import spatial_correlation_matrix
+from q01 import array_manifold_vector_linear
+
+
+def compute_spatial_spectrum(z, p, fs):
+    F, T, Z = stft(z, fs, window="hann", nperseg=1024, noverlap=512)
+    angles = np.linspace(0, 360, 361)
+    P = np.zeros((31, len(angles)))
+    R_z = spatial_correlation_matrix(Z)
+    tau = np.array([0, 10 / fs, 20 / fs])
+
+    for i, theta in enumerate(angles):
+        w = 1 / 3 * np.exp(-1j * 2 * np.pi * theta * tau)
+        for f_bin in range(20, 31):
+            P[f_bin, i] = np.abs(np.conj(w.T).dot(R_z[f_bin]).dot(w))
+
+    plt.figure(figsize=(10, 6))
+    for f_bin in range(20, 31):
+        plt.plot(angles, 20 * np.log10(P[f_bin, :]), label=f"Frequency bin {f_bin}")
+    plt.xlabel("Angle [degrees]")
+    plt.ylabel("20log10|P(θ)|")
+    plt.title("Spatial Spectrum")
+    plt.legend()
+    plt.show()
+
+
+fs = 16000
+T = 1
+N = fs * T
+t = np.linspace(0, T, N, endpoint=False)
+f = 440
+a = 1
+
+s = a * np.sin(2 * np.pi * f * t)
+
+np.random.seed(0)
+noise = np.random.normal(0, 1, N)
+SNR = 10
+noise_amplitude = a / (10 ** (SNR / 20))
+epsilon = noise_amplitude * noise
+
+x1 = s + epsilon
+x2 = np.roll(s, 10) + epsilon
+x3 = np.roll(s, 20) + epsilon
+x = np.array([x1, x2, x3])
+
+
+d = 0.05
+M = 3
+linear_coords = [(m * d, 0) for m in range(M)]
+
+compute_spatial_spectrum(x, linear_coords, fs)