diff --git a/yyamamoto/chapter05/q01.py b/yyamamoto/chapter05/q01.py new file mode 100644 index 0000000..3b9cae1 --- /dev/null +++ b/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] \ No newline at end of file diff --git a/yyamamoto/chapter05/q02.py b/yyamamoto/chapter05/q02.py new file mode 100644 index 0000000..d669097 --- /dev/null +++ b/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] \ No newline at end of file diff --git a/yyamamoto/chapter05/q03.py b/yyamamoto/chapter05/q03.py new file mode 100644 index 0000000..9577919 --- /dev/null +++ b/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] +# 合ってる! \ No newline at end of file diff --git a/yyamamoto/chapter05/q04.py b/yyamamoto/chapter05/q04.py new file mode 100644 index 0000000..7f2df2d --- /dev/null +++ b/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]]] \ No newline at end of file diff --git a/yyamamoto/chapter05/q05.py b/yyamamoto/chapter05/q05.py new file mode 100644 index 0000000..4845786 --- /dev/null +++ b/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]] \ No newline at end of file diff --git a/yyamamoto/chapter05/q06.py b/yyamamoto/chapter05/q06.py new file mode 100644 index 0000000..89378c6 --- /dev/null +++ b/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") + + + diff --git a/yyamamoto/chapter05/q06_graph.png b/yyamamoto/chapter05/q06_graph.png new file mode 100644 index 0000000..1211022 Binary files /dev/null and b/yyamamoto/chapter05/q06_graph.png differ diff --git a/yyamamoto/chapter05/q07.py b/yyamamoto/chapter05/q07.py new file mode 100644 index 0000000..40dc945 --- /dev/null +++ b/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') + diff --git a/yyamamoto/chapter05/q07_graph.png b/yyamamoto/chapter05/q07_graph.png new file mode 100644 index 0000000..19df6d5 Binary files /dev/null and b/yyamamoto/chapter05/q07_graph.png differ diff --git a/yyamamoto/chapter05/q08.py b/yyamamoto/chapter05/q08.py new file mode 100644 index 0000000..b6142be --- /dev/null +++ b/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) \ No newline at end of file diff --git a/yyamamoto/chapter05/q08_graph.png b/yyamamoto/chapter05/q08_graph.png new file mode 100644 index 0000000..eabc5a5 Binary files /dev/null and b/yyamamoto/chapter05/q08_graph.png differ diff --git a/yyamamoto/chapter05/q09.py b/yyamamoto/chapter05/q09.py new file mode 100644 index 0000000..65830a9 --- /dev/null +++ b/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') \ No newline at end of file diff --git a/yyamamoto/chapter05/q09_graph.png b/yyamamoto/chapter05/q09_graph.png new file mode 100644 index 0000000..25d0b5a Binary files /dev/null and b/yyamamoto/chapter05/q09_graph.png differ diff --git a/yyamamoto/chapter05/q10.py b/yyamamoto/chapter05/q10.py new file mode 100644 index 0000000..8aa335b --- /dev/null +++ b/yyamamoto/chapter05/q10.py @@ -0,0 +1,17 @@ +import numpy as np +import matplotlib.pyplot as plt +from q05 import stft +from q04 import correlation_matrix + +def spatial_spectrum(z): + M = z.shape[0] + L = 1024 + S = 512 + window = np.hanning(L) + T, F = stft(L, S, window, z[0]).shape + Z = np.zeros((M, F, T), dtype=complex) + for m in range(M): + Z[m] = stft(L, S, w, z).T + corr_matrix = correlation_matrix(Z) + +