add chapter05後半

eimamura/chapter05/q06.py

@@ -0,0 +1,37 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def calculate_amplitude_from_snr(s, snr):
+    a = np.linalg.norm(s)
+    x = np.random.randn(round(len(s)))
+    x = a * x / (10 ** (snr / 20) * np.linalg.norm(x))
+    return x
+
+
+if __name__ == "__main__":
+    snr = 10
+    s = 1
+    fs = 16000
+    a = 1
+    f = 440
+
+    t = np.arange(0, s, 1 / fs)
+    y = a * np.sin(2 * np.pi * f * t)
+
+    s1 = a * np.sin(2 * np.pi * f * t)
+    s2 = np.hstack((np.zeros(10), s1[:-10]))
+    s3 = np.hstack((np.zeros(20), s1[:-20]))
+
+    white_noise = calculate_amplitude_from_snr(y, snr)
+
+    x1 = s1 + white_noise
+    x2 = s2 + white_noise
+    x3 = s3 + white_noise
+
+    plt.plot(x1, label="x1")
+    plt.plot(x2, label="x2")
+    plt.plot(x3, label="x3")
+    plt.xlim(0, 0.01 * fs)
+    plt.legend(["x1[n]", "x2[n]", "x3[n]"])
+    plt.show()

eimamura/chapter05/q07.py

@@ -0,0 +1,118 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def pad(x, L, S):
+    x_pad = np.pad(x, [L - S, L - S])
+    re = np.mod(x_pad.size, S)
+    if re != 0:
+        x_pad = np.pad(x_pad, [0, S - re])
+    return x_pad
+
+
+def frame_div(x, L, S):
+    x_pad = pad(x, L, S)
+    T = int(np.floor((x_pad.size - L) / S)) + 1
+    x_t = np.array([x_pad[t * S : t * S + L] for t in range(T)])
+    return x_t
+
+
+def stft(x, L, S, wnd):
+    x_t = frame_div(x, L, S)
+    T = len(x_t)
+    X = np.array([np.fft.rfft(x_t[t] * wnd) for t in range(T)], dtype="complex")
+    return X.T
+
+
+def composite(S, w):
+    L = len(w)
+    Q = int(L / S)
+    ws = np.zeros(L)
+    a = 0
+    for l in range(L):
+        for m in range(-(Q - 1), Q):
+            if (l - m * S) >= 0 and L > (l - m * S):
+                a += w[l - m * S] ** 2
+        ws[l] = w[l] / a
+        a = 0
+
+    return ws
+
+
+def istft(S, X):
+    F = X.shape[0]
+    T = X.shape[1]
+    N = 2 * (F - 1)
+    M = S * (T - 1) + N
+
+    com = np.hamming(N)
+    com1 = composite(S, com)
+    x = np.zeros(M)
+    z = np.zeros((T, N))
+    for t in range(T):
+        for n in range(N):
+            z[t, n] = np.fft.irfft(X[:, t])[n]
+            x[t * S + n] += com1[n] * z[t, n]
+
+    return x
+
+
+def calculate_amplitude_from_snr(s, snr):
+    a = np.linalg.norm(s)
+    x = np.random.randn(round(len(s)))
+    x = a * x / (10 ** (snr / 20) * np.linalg.norm(x))
+    return x
+
+
+if __name__ == "__main__":
+    snr = 10
+    s = 1
+    fs = 16000
+    a = 1
+    f = 440
+    L = 1024
+    S = 512
+
+    t = np.arange(0, s, 1 / fs)
+    y = a * np.sin(2 * np.pi * f * t)
+
+    s1 = a * np.sin(2 * np.pi * f * t)
+    s2 = np.hstack((np.zeros(10), s1[:-10]))
+    s3 = np.hstack((np.zeros(20), s1[:-20]))
+    white_noise = calculate_amplitude_from_snr(y, snr)
+    x1 = s1 + white_noise
+    x2 = s2 + white_noise
+    x3 = s3 + white_noise
+
+    wnd = np.hamming(L)
+    X1 = stft(x1, L, S, wnd)
+    X2 = stft(x2, L, S, wnd)
+    X3 = stft(x3, L, S, wnd)
+
+    X = np.stack([X1, X2, X3])
+    M, F, T = X.shape
+
+    f = (fs / 2) / (F - 1) * np.arange(F)
+    tau1 = 0
+    tau2 = 10 / fs
+    tau3 = 20 / fs
+    w = (
+        1
+        / 3
+        * np.array(
+            [
+                np.exp(-1j * 2 * np.pi * f * tau1),
+                np.exp(-1j * 2 * np.pi * f * tau2),
+                np.exp(-1j * 2 * np.pi * f * tau3),
+            ]
+        )
+    )
+
+    Y = np.sum(np.conjugate(w[:, :, None]) * X, axis=0)
+
+    a = istft(S, Y)
+    plt.plot(x1)
+    plt.plot(a[L - S :])
+    plt.xlim([0, 0.01 * fs])
+    plt.legend(["obs,(x1)", "enh"])
+    plt.show()

eimamura/chapter05/q08.py

@@ -0,0 +1,61 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def general_array_manifold_vector(p, theta, fs):
+    c = 334
+    u = np.array([np.sin(np.radians(theta)), np.cos(np.radians(theta)), 0])
+    M = len(p)
+    a = np.zeros(M, dtype="complex")
+
+    for m in range(M):
+        a[m] = np.exp(1j * 2 * np.pi * fs / c * u @ p[m])
+
+    return a
+
+
+def beam_pattern(w, p, fs):
+    """ビームパターン
+    Args:
+        w(array-like): 周波数領域におけるビームフォーマのフィルタ
+        p(array-like): マイクアレイの座標
+        fs(int): サンプリング周波数
+    Return:
+        ビームパターンを描画
+    """
+    F = w.shape[0]
+
+    deg = np.arange(360)
+    f = np.arange(F) * (fs / 2) / (F - 1)
+    psi = np.zeros((F, 360), dtype="complex")
+
+    for _f in range(F):
+        for theta in range(360):
+            a = general_array_manifold_vector(p, theta, f[_f])
+            psi[_f, theta] = np.conjugate(w[_f]) @ a
+
+    plt.pcolormesh(deg, f, 20 * np.log10(abs(psi)))
+    plt.colorbar()
+    plt.show()
+
+    return
+
+
+if __name__ == "__main__":
+    M = 3
+    d = 0.05
+    F = 512
+    fs = 16000
+    p_linear = []
+    for m in range(M):
+        p_linear.append([((m - 1) - (M - 1) / 2) * d, 0, 0])
+    p_linear = np.array(p_linear)
+
+    beam_angle = 0
+    f = (fs / 2) / (F - 1) * np.arange(F)
+    w = []
+    for _f in f:
+        w.append(general_array_manifold_vector(p_linear, beam_angle, _f))
+    w = np.array(w) / M
+
+    beam_pattern = beam_pattern(w, p_linear, fs)

eimamura/chapter05/q09.py

@@ -0,0 +1,56 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def general_array_manifold_vector(p, theta, fs):
+    c = 334
+    u = np.array([np.sin(np.radians(theta)), np.cos(np.radians(theta)), 0])
+    M = len(p)
+    a = np.zeros(M, dtype="complex")
+
+    for m in range(M):
+        a[m] = np.exp(1j * 2 * np.pi * fs / c * u @ p[m])
+
+    return a
+
+
+def beam_pattern(w, p, fs):
+    F = w.shape[0]
+
+    deg = np.arange(360)
+    f = np.arange(F) * (fs / 2) / (F - 1)
+    psi = np.zeros((F, 360), dtype="complex")
+
+    for _f in range(F):
+        for theta in range(360):
+            a = general_array_manifold_vector(p, theta, f[_f])
+            psi[_f, theta] = np.conjugate(w[_f]) @ a
+
+    plt.pcolormesh(deg, f, 20 * np.log10(abs(psi)))
+    plt.colorbar()
+    plt.show()
+
+    return
+
+
+if __name__ == "__main__":
+    M = 3
+    D = np.array([0.02, 0.05, 0.10])
+    F = 1000
+    fs = 16000
+
+    for d in D:
+        w = np.zeros((F, M), dtype="complex")
+        p_linear = []
+        for m in range(M):
+            p_linear.append([((m - 1) - (M - 1) / 2) * d, 0, 0])
+        p_linear = np.array(p_linear)
+
+        beam_angle = 0
+        f = (fs / 2) / (F - 1) * np.arange(F)
+        w = []
+        for _f in f:
+            w.append(general_array_manifold_vector(p_linear, beam_angle, _f))
+        w = np.array(w) / M
+
+        beam_pattern(w, p_linear, fs)

eimamura/chapter05/q10.py

@@ -0,0 +1,106 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def general_array_manifold_vector(p, theta, fs):
+    c = 334
+    u = np.array([np.sin(np.radians(theta)), np.cos(np.radians(theta)), 0])
+    M = len(p)
+    a = np.zeros(M, dtype="complex")
+
+    for m in range(M):
+        a[m] = np.exp(1j * 2 * np.pi * fs / c * u @ p[m])
+
+    return a
+
+
+def spatial_correlation_matrix(X):
+    M, F, T = np.shape(X)
+    xft = np.zeros((F, M, T), dtype=np.complex64)
+    for f in range(F):
+        for m in range(M):
+            xft[f][m] = X[m][f]
+
+    R = np.zeros((F, M, M), dtype=np.complex64)
+    for f in range(F):
+        R[f] = np.dot(xft[f], np.conjugate(xft[f].T)) / T
+    return R
+
+
+def pad(x, L, S):
+    x_pad = np.pad(x, [L - S, L - S])
+    re = np.mod(x_pad.size, S)
+    if re != 0:
+        x_pad = np.pad(x_pad, [0, S - re])
+    return x_pad
+
+
+def frame_div(x, L, S):
+    x_pad = pad(x, L, S)
+    T = int(np.floor((x_pad.size - L) / S)) + 1
+    x_t = np.array([x_pad[t * S : t * S + L] for t in range(T)])
+    return x_t
+
+
+def stft(x, L, S, wnd):
+    x_t = frame_div(x, L, S)
+    T = len(x_t)
+    X = np.array([np.fft.rfft(x_t[t] * wnd) for t in range(T)], dtype="complex")
+    return X.T
+
+
+def calculate_amplitude_from_snr(s, snr):
+    a = np.linalg.norm(s)
+    x = np.random.randn(round(len(s)))
+    x = a * x / (10 ** (snr / 20) * np.linalg.norm(x))
+    return x
+
+
+def spatial_spectrum(z, f):
+    L = 1024
+    S = 512
+    d = 0.05
+    fs = 16000
+    p = np.array([[-d, 0, 0], [0, 0, 0], [d, 0, 0]])
+    win = np.hanning(L)
+    M = z.shape[0]
+
+    Z = []
+    for m in range(M):
+        Z.append(stft(z[m], L, S, win))
+    Z = np.array(Z)
+
+    R = spatial_correlation_matrix(Z)
+    thetas = np.arange(360)
+    w = []
+    for theta in thetas:
+        a = general_array_manifold_vector(p, theta, fs / 2 / (Z.shape[1] - 1) * f)
+        w.append(a / (a.conj() @ a))
+    w = np.array(w)
+
+    P = []
+    for theta in thetas:
+        P.append(np.dot(w[theta].conj(), R[f]) @ w[theta])
+    return np.array(P)
+
+
+if __name__ == "__main__":
+    snr = 10
+    s = 1
+    fs = 16000
+    a = 1
+    f = 440
+
+    t = np.arange(0, s, 1 / fs)
+    y = a * np.sin(2 * np.pi * f * t)
+
+    white_noise = calculate_amplitude_from_snr(y, snr)
+
+    x = np.zeros((3, len(y)))
+    for i in range(len(x)):
+        x[i] = np.roll(y, i * 10) + white_noise
+
+    for i in range(20, 30):
+        plt.plot(np.arange(360), 20 * np.log10(np.abs(spatial_spectrum(x, i))))
+        plt.title("bin : " + str(i))
+        plt.show()