5章前半を追加

yyamamoto/chapter05/q01.py

@@ -0,0 +1,21 @@
+import numpy as np
+
+def straight_array_manifold_vector(d, M, theta, f):
+    c = 334
+    u = np.array([np.sin(theta), np.cos(theta), 0]).T
+    p = np.zeros((M, 3))
+    a = np.zeros(M, dtype=complex)
+    for m in range(1, M + 1):
+        p[m - 1] = 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.T, p[m - 1]))
+    return a
+
+# 確認
+d = 0.05
+M = 3
+theta = np.pi / 4
+f = 1000
+print(straight_array_manifold_vector(d, M, theta, f))
+
+# 結果
+# [0.7868537-0.61713958j 1.       +0.j         0.7868537+0.61713958j]

yyamamoto/chapter05/q02.py

@@ -0,0 +1,21 @@
+import numpy as np
+
+def circular_array_manifold_vector(r, M, theta, f):
+    c = 334
+    u = np.array([np.sin(theta), np.cos(theta), 0]).T
+    p = np.zeros((M, 3))
+    a = np.zeros(M, dtype=complex)
+    for m in range(1, M + 1):
+        p[m - 1] = 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.T, p[m - 1]))
+    return a
+
+# 確認
+d = 0.05
+M = 3
+theta = np.pi / 4
+f = 1000
+print(circular_array_manifold_vector(d, M, theta, f))
+
+# 結果
+# [0.7868537 +0.61713958j 0.97051349+0.2410468j  0.6148926 -0.78861086j]

yyamamoto/chapter05/q03.py

@@ -0,0 +1,35 @@
+import numpy as np
+
+def array_manifold_vector(p, M, theta, f):
+    M = p.shape[0]
+    c = 334
+    u = np.array([np.sin(theta), np.cos(theta), 0]).T
+    a = np.zeros(M, dtype=complex)
+    for m in range(1, M + 1):
+        a[m - 1] = np.exp(1j * 2 * np.pi * f / c * np.dot(u.T, p[m - 1]))
+    return a
+
+# 1の確認
+d = 0.05
+M = 3
+theta = np.pi / 4
+f = 1000
+p = np.zeros((M, 3))
+for m in range(1, M + 1):
+    p[m - 1] = np.array([((m - 1) - (M - 1) / 2) * d, 0, 0]).T
+print(array_manifold_vector(p, M, theta, f))
+
+# 1の結果
+# [0.7868537-0.61713958j 1.       +0.j         0.7868537+0.61713958j]
+# 合ってる！
+
+
+# 2の確認
+r = 0.05
+for m in range(1, M + 1):
+    p[m - 1] = np.array([r * np.sin(2 * np.pi / M * (m - 1)), r * np.cos(2 * np.pi / M * (m - 1)), 0]).T
+print(array_manifold_vector(p, M, theta, f))
+
+# 2の結果
+# [0.7868537 +0.61713958j 0.97051349+0.2410468j  0.6148926 -0.78861086j]
+# 合ってる！

yyamamoto/chapter05/q04.py

@@ -0,0 +1,31 @@
+import numpy as np
+
+def 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 = np.array([X[:, f, t]]).T
+            R[f] += np.dot(x_ft, np.conjugate(x_ft.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])
+
+print(correlation_matrix(X))
+
+
+# 結果の確認
+# [[[ 4.+0.j  5.+0.j]
+#   [ 5.+0.j 21.+0.j]]
+
+#  [[16.+0.j  6.+0.j]
+#   [ 6.+0.j  5.+0.j]]
+
+#  [[36.+0.j  3.+0.j]
+#   [ 3.+0.j  3.+0.j]]]

yyamamoto/chapter05/q05.py

@@ -0,0 +1,57 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q04 import correlation_matrix
+
+def zero_padding(L, S, x):
+    N = len(x)
+    length = L - S + N
+    if length % S != 0:
+        length += S - (length % S)
+    length += L - S
+    ans = np.zeros(length)
+    ans[L-S:L-S+N] = x
+
+    return ans
+
+def split_frame(L, S, x):
+    x_tilde = zero_padding(L, S, x)
+    N = len(x_tilde)
+    T = (N - L) // S + 1
+
+    ans = np.zeros((T, L))
+    for t in range(T):
+        ans[t] = x_tilde[t*S:t*S+L]
+    return ans
+
+def stft(L, S, w, x):
+    x_splited = split_frame(L, S, x)
+    T = x_splited.shape[0]
+    x_stft = np.zeros((T, L//2+1), dtype=complex)
+    for t in range(T):
+        x_splited[t] = x_splited[t] * w
+    for t in range(T):
+        x_stft[t] = np.fft.rfft(x_splited[t])
+    return x_stft
+
+
+
+fs = 16000  # サンプリング周波数 Hz
+sec = 5       # 信号長 s
+
+t = np.arange(0, sec, 1/fs) # サンプリング点の配列
+white_noise1 = 2 * (np.random.normal(size = fs * sec)) - 1
+white_noise2 = 2 * (np.random.normal(size = fs * sec)) - 1
+
+
+L = 512
+S = 256
+w = np.hanning(L)
+
+X1 = stft(L, S, w, white_noise1)
+X2 = stft(L, S, w, white_noise2)
+X = np.stack([X1, X2])
+print(correlation_matrix(X)[100].real)
+
+# 結果の確認
+# [[272642.81798742  60046.86665743]
+#  [ 60046.86665743 228203.33404359]]

yyamamoto/chapter05/q06.py

@@ -0,0 +1,51 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def white_noise(s, sn):
+    N = len(s)
+    x = 2 * (np.random.rand(N)) - 1     # ホワイトノイズの作成
+
+    sum_s = 0
+    sum_x = 0
+    for n in range(N):
+        sum_s += s[n] * s[n]
+        sum_x += x[n] * x[n]
+
+    mul = pow(np.exp(-sn/10) * sum_s / sum_x, 0.5)
+    x = x * mul
+
+    # return x
+    ans = s + x
+    return ans
+
+fs = 16000
+time = 1
+f = 440
+t = np.arange(-20, fs * time) / fs      # 20サンプル分手前に作っておく
+s = np.sin(2 * np.pi * f * t)
+epsilon = white_noise(s[20::], 10)
+x1 = s[20:] + epsilon
+x2 = s[10:len(s) - 10] + epsilon
+x3 = s[:len(s) - 20] + epsilon
+x = np.array([x1, x2, x3])
+
+fig = plt.figure()
+plt.subplot(3, 1, 1)
+plt.xlim([0.5, 0.51])
+plt.plot(t[20:], x1)
+plt.xlabel('time[s]')
+plt.title("x1")
+plt.subplot(3, 1, 2)
+plt.xlim([0.5, 0.51])
+plt.plot(t[10:len(s) - 10], x2)
+plt.xlabel('time[s]')
+plt.title("x2")
+plt.subplot(3, 1, 3)
+plt.xlim([0.5, 0.51])
+plt.plot(t[:len(s) - 20], x3)
+plt.xlabel('time[s]')
+plt.title("x3")
+plt.savefig("./yyamamoto/chapter05/q06_graph.png")
+
+
+

yyamamoto/chapter05/q07.py

@@ -0,0 +1,75 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q05 import stft
+from q06 import white_noise
+
+def window(S, w):
+    L = len(w)
+    Q = L // S
+    w_s = np.zeros(L)
+    for l in range(L):
+        den = 0
+        for m in range(-(Q-1), Q):
+            if 0 <= l - m * S & l - m * S < L:
+                den += w[l-m*S] ** 2
+        w_s[l] = w[l] / den
+    return w_s
+
+def istft(S, X):
+    # 手順1
+    F, T = X.shape
+    N = 2 * (F - 1)
+    M = S * (T - 1) + N
+    # 手順2
+    x = np.zeros(M)
+    # 手順3
+    z = np.zeros((2*(F-1), T))
+    for t in range(T):
+        z[:,t] = np.fft.irfft(X[:,t])
+    # 手順4
+    w = window(S, np.hanning(N))
+
+    n = np.arange(N)    # ファンシーインデックス
+    for t in range(T):
+        x[t*S+n] = x[t*S+n] + w[n] * z[n,t]
+    return x
+
+fs = 16000
+time = 1
+f = 440
+t = np.arange(-20, fs * time) / fs      # 20サンプル分手前に作っておく
+s = np.sin(2 * np.pi * f * t)
+epsilon = white_noise(s[20::], 10)
+x1 = s[20:] + epsilon
+x2 = s[10:len(s) - 10] + epsilon
+x3 = s[:len(s) - 20] + epsilon
+
+# ここからq07
+L = 1024
+S = 512
+han = np.hanning(L)
+X1 = stft(L, S, han, x1).T
+X2 = stft(L, S, han, x2).T
+X3 = stft(L, S, han, x3).T
+F, T = X1.shape
+Y = np.zeros((X1.shape), dtype=complex)
+for f_idx in range(F):
+    f = fs / 2 / (F - 1) * f_idx
+    w_f = np.exp(np.array([0, -1j * 2 * np.pi * f * 10 / fs, -1j * 2 * np.pi * f * 20 / fs])) / 3
+    for t in range(T):
+        x_ft = np.array([[X1[f_idx][t], X2[f_idx][t], X3[f_idx][t]]]).T
+        Y[f_idx][t] = (np.conjugate(w_f) @ x_ft)[0]
+
+print(Y)
+Y_istft = istft(S, Y)
+
+print("\n\n\n\n")
+print(Y_istft)
+
+t = np.arange(len(Y_istft)) / fs
+fig = plt.figure()
+plt.plot(t, Y_istft)
+plt.xlim([0.5, 0.51])
+plt.xlabel("time[s]")
+fig.savefig('./yyamamoto/chapter05/q07_graph.png')
+

yyamamoto/chapter05/q08.py

@@ -0,0 +1,42 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import array_manifold_vector
+
+def show_beam_pattern(w, p, fs):
+    # 手順a
+    F = w.shape[0]
+    M = p.shape[0]
+    a = np.zeros((F, 361, M), dtype=complex)
+    for f_idx in range(F):
+        f = fs / 2 / (F - 1) * f_idx
+        for theta_idx in range(361):
+            theta = theta_idx / 180 * np.pi     # ラジアンに変換
+            a[f_idx][theta_idx] = array_manifold_vector(p, M, theta, f)
+    # 手順b
+    psi = np.zeros((F, 361), dtype=complex)
+    for f_idx in range(F):
+        for theta_idx in range(361):
+            psi[f_idx][theta_idx] = w[f_idx] @ a[f_idx][theta_idx].T
+
+    psi_db = 20 * np.log10(np.abs(psi)).T
+    thetas = np.arange(361)
+    freqs = np.arange(F) * fs / 2 / (F - 1)
+    fig = plt.figure()
+    plt.pcolormesh(thetas, freqs, psi_db.T)
+    plt.xlabel('angle[deg]')
+    plt.ylabel('frequency[Hz]')
+    fig.savefig('./yyamamoto/chapter05/q08_graph.png')
+
+d = 0.05
+M = 3
+F = 1000
+fs = 16000
+p = np.zeros((M, 3))
+w = np.zeros((F, 3), dtype=complex)
+for m in range(1, M + 1):
+    p[m - 1] = np.array([((m - 1) - (M - 1) / 2) * d, 0, 0]).T
+tau = np.array([0, 10 / fs, 20 / fs])
+for f in range(F):
+    w[f] = np.exp(-1j * 2 * np.pi * f * tau) / 3
+
+show_beam_pattern(w, p, fs)

yyamamoto/chapter05/q09.py

@@ -0,0 +1,49 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import array_manifold_vector
+
+def show_beam_pattern(w, p, fs):
+    # 手順a
+    F = w.shape[0]
+    M = p.shape[0]
+    a = np.zeros((F, 361, M), dtype=complex)
+    for f_idx in range(F):
+        f = fs / 2 / (F - 1) * f_idx
+        for theta_idx in range(361):
+            theta = theta_idx / 180 * np.pi     # ラジアンに変換
+            a[f_idx][theta_idx] = array_manifold_vector(p, M, theta, f)
+    # 手順b
+    psi = np.zeros((F, 361), dtype=complex)
+    for f_idx in range(F):
+        for theta_idx in range(361):
+            psi[f_idx][theta_idx] = w[f_idx] @ a[f_idx][theta_idx].T
+
+    psi_db = 20 * np.log10(np.abs(psi)).T
+    thetas = np.arange(361)
+    freqs = np.arange(F) * fs / 2 / (F - 1)
+
+    plt.pcolormesh(thetas, freqs, psi_db.T)
+    plt.xlabel('angle[deg]')
+    plt.ylabel('frequency[Hz]')
+
+
+
+
+M = 3
+F = 1000
+fs = 16000
+p = np.zeros((M, 3))
+w = np.zeros((F, 3), dtype=complex)
+ds = [0.02, 0.05, 0.1]
+fig = plt.figure()
+for i in range(1, 4):
+    d = ds[i - 1]
+    for m in range(1, M + 1):
+        p[m - 1] = np.array([((m - 1) - (M - 1) / 2) * d, 0, 0]).T
+    tau = np.array([0, 10 / fs, 20 / fs])
+    for f in range(F):
+        w[f] = np.exp(-1j * 2 * np.pi * f * tau) / 3
+    plt.subplot(1, 3, i)
+    show_beam_pattern(w, p, fs)
+
+fig.savefig('./yyamamoto/chapter05/q09_graph.png')