diff --git a/dsugawara/chapter05/q01.py b/dsugawara/chapter05/q01.py new file mode 100644 index 0000000..ce6bc22 --- /dev/null +++ b/dsugawara/chapter05/q01.py @@ -0,0 +1,34 @@ +import numpy as np + + +def linear_array_vector(d, M, theta, f): + """直線状アレイのアレイマニフォールドベクトルを求める + Args: + d (double): アレイ間隔 + M (int): マイク数 + theta (int): 音源方向（y軸から反時計回りが正の向き） + f (int): 周波数 + Return: + amv: アレイマニフォールドベクトル + """ + c = 334.0 + theta = np.pi / 2 + 2 * np.pi * (theta / 360) + u = np.array([np.sin(theta), np.cos(theta), 0]).T + + p = np.zeros([M, 3]) + for m in range(1, M + 1): + idx_m = m - 1 + p[idx_m, 0] = ((m - 1) - (M - 1) / 2) * d + p = p.T + amv = np.empty([u.T.shape[0], p.shape[1]], dtype="complex") + amv = np.exp(2j * np.pi * f / c * np.dot(u.T, p)) + + return amv + + +if __name__ == "__main__": + d = 0.05 + M = 3 + theta = 45 + f = 1000 + print(linear_array_vector(d, M, theta, f)) diff --git a/dsugawara/chapter05/q02.py b/dsugawara/chapter05/q02.py new file mode 100644 index 0000000..38b8481 --- /dev/null +++ b/dsugawara/chapter05/q02.py @@ -0,0 +1,35 @@ +import numpy as np + + +def circular_array_vector(r, M, theta, f): + """円状アレイのアレイマニフォールドベクトルを求める + Args: + r (double): アレイ半径 + M (int): マイク数 + theta (int): 音源方向（y軸から反時計回りが正の向き） + f (int): 周波数 + Return: + amv: アレイマニフォールドベクトル + """ + c = 334.0 + theta = np.pi / 2 + 2 * np.pi * (theta / 360) + u = np.array([np.sin(theta), np.cos(theta), 0]).T + + p = np.zeros([M, 3]) + for m in range(1, M + 1): + idx_m = m - 1 + p[idx_m, 0] = r * np.sin(2 * np.pi / M * (m - 1)) + p[idx_m, 1] = r * np.cos(2 * np.pi / M * (m - 1)) + p = p.T + amv = np.empty([u.T.shape[0], p.shape[1]], dtype="complex") + amv = np.exp(2j * np.pi * f / c * np.dot(u.T, p)) + + return amv + + +if __name__ == "__main__": + r = 0.05 + M = 3 + theta = 45 + f = 1000 + print(circular_array_vector(r, M, theta, f)) diff --git a/dsugawara/chapter05/q03.py b/dsugawara/chapter05/q03.py new file mode 100644 index 0000000..d369bc0 --- /dev/null +++ b/dsugawara/chapter05/q03.py @@ -0,0 +1,54 @@ +import numpy as np + + +def array_vector(crd, theta, f): + """アレイの座標からアレイマニフォールドベクトルを求める + Args: + crd (ndarray): アレイの座標 + theta (int): 音源方向（y軸から反時計回りが正の向き） + f (int): 周波数 + Return: + amv: アレイマニフォールドベクトル + """ + M = crd.shape[0] + c = 334.0 + theta = np.pi / 2 + 2 * np.pi * (theta / 360) + u = np.array([np.sin(theta), np.cos(theta), 0]).T + + p = np.zeros([M, 3]) + for m in range(1, M + 1): + idx_m = m - 1 + r = crd[idx_m, 0] ** 2 + crd[idx_m, 1] ** 2 + if r == 0: + sin, cos = 0, 0 + else: + sin = crd[idx_m, 1] / r + cos = crd[idx_m, 0] / r + p[idx_m, 0] = r * sin + p[idx_m, 1] = r * cos + p = p.T + amv = np.empty([u.T.shape[0], p.shape[1]], dtype="complex") # a[1, M]になるはず + amv = np.exp(2j * np.pi * f / c * np.dot(u.T, p)) + + return amv + + +if __name__ == "__main__": + d = 0.05 + theta = 45 + f = 1000 + + # 直線状アレイの確認 + liner = np.array([[-d, 0, 0], [0, 0, 0], [d, 0, 0]]) + print(array_vector(liner, theta, f)) + + # 円状アレイの確認 + r = 0.05 + circular = np.array( + [ + [r, 0, 0], + [r * np.cos(2 * np.pi / 3), r * np.sin(2 * np.pi / 3), 0], + [r * np.cos(2 * np.pi * 2 / 3), r * np.sin(2 * np.pi * 2 / 3), 0], + ] + ) + print(array_vector(circular, theta, f)) diff --git a/dsugawara/chapter05/q04.py b/dsugawara/chapter05/q04.py new file mode 100644 index 0000000..0c1f1f1 --- /dev/null +++ b/dsugawara/chapter05/q04.py @@ -0,0 +1,28 @@ +import numpy as np + + +def scm(X): + """空間相関行列（spatial correation matrix）を求める + Args: + X (ndarray): 複素数行列 + Return: + R: 空間相関行列 + """ + M, F, T = X.shape + + R = np.empty([F, M, M], dtype="complex") + x = np.empty([F, M, T], dtype="complex") + for f in range(0, F): + for m in range(0, M): + x[f, m] = X[m, f] + for f in range(0, F): + R[f] = np.dot(x[f], np.conjugate(x[f].T)) / T + return R + + +if __name__ == "__main__": + 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.array([X1, X2]) + + print(scm(X)) diff --git a/dsugawara/chapter05/q05.py b/dsugawara/chapter05/q05.py new file mode 100644 index 0000000..d792b3a --- /dev/null +++ b/dsugawara/chapter05/q05.py @@ -0,0 +1,103 @@ +import numpy as np + + +def padding(L, S, x): + """N点の信号xにゼロ詰めを行う + Args: + L(int): 窓幅 + S(int): シフト幅 + x(ndarray): 入力信号 + Return + x_pad: ゼロ詰め後の信号 + """ + + x_pad = np.pad(x, (L - S, L - S)) + temp = S - len(x_pad) % S + x_pad = np.pad(x_pad, (0, temp)) + return x_pad + + +def flame_div(L, S, x): + """N点の信号xをフレーム分割 + Args: + L(int): 窓幅 + S(int): シフト幅 + x(ndarray): 入力信号 + Return + x_div: フレーム分割したx + """ + x_pad = padding(L, S, x) + N = len(x_pad) + T = ((N - L) // S) + 1 + x_div = np.empty([T, L]) + for t in range(0, T): + x_div[t] = x_pad[t * S : t * S + L] + + return x_div + + +def my_stft(L, S, x, win): + """短時間フーリエ変換 + Args: + L(int): 窓幅 + S(int): シフト幅 + x(ndarray): 入力信号 + w(ndarray): 窓関数 + Return + X: 変換後の信号 + """ + x_div = flame_div(L, S, x) + T = x_div.shape[0] + X = np.empty([L // 2 + 1, T], dtype="complex") + for t in range(0, T): + X[:, t] = np.fft.rfft(x_div[t, :] * win) + + return X + + +def scm(X): + """空間相関行列（spatial correation matrix）を求める + Args: + X (ndarray): 複素数行列 + Return: + R: 空間相関行列 + """ + M, F, T = X.shape + + R = np.empty([F, M, M], dtype="complex") + x = np.empty([F, M, T], dtype="complex") + for f in range(0, F): + for m in range(0, M): + x[f, m] = X[m, f] + for f in range(0, F): + R[f] = np.dot(x[f], np.conjugate(x[f].T)) / T + return R + + +if __name__ == "__main__": + 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.array([X1, X2]) + + print(scm(X)) + + +if __name__ == "__main__": + + A = 1.0 + fs = 16000 + sec = 5 + + nA = A * np.random.randn(round(fs * sec)) + nB = A * np.random.rand(round(fs * sec)) + + L = 512 # 窓長 + S = 256 # シフト幅 + win = np.hanning(L) # 窓関数 + + nA_stft = my_stft(L, S, nA, win) + nB_stft = my_stft(L, S, nB, win) + + R = scm(np.array([nA_stft, nB_stft])) + + print(R[100].real) diff --git a/dsugawara/chapter05/q06.py b/dsugawara/chapter05/q06.py new file mode 100644 index 0000000..3ae8a2c --- /dev/null +++ b/dsugawara/chapter05/q06.py @@ -0,0 +1,44 @@ +import numpy as np +import matplotlib.pyplot as plt + + +def make_noise(signal, is_SNR): + """SNRを調整した雑音を生成する + Args: + signal(ndarray): 信号 + is_SNR(double): SNR + Return: + ndarray: 雑音 + """ + noise = np.random.randn(len(signal)) + + noise = noise / np.sqrt(np.sum(noise**2)) + noise = noise * np.sqrt(np.sum(signal**2)) + noise = noise * 10 ** (-1 * is_SNR / 20) + + return noise + + +if __name__ == "__main__": + A = 1.0 + fs = 16000 + sec = 1 + f = 440 + t = np.arange(0, fs * sec) / fs + + s = A * np.sin(2 * np.pi * f * t) + + noise = make_noise(s, 10) + + x1 = s + noise + s2 = np.pad(s, (10, 0)) + x2 = s2[0 : len(s2) - 10] + noise + s3 = np.pad(s, (20, 0)) + x3 = s3[0 : len(s3) - 20] + noise + + plt.figure(figsize=[6.0, 4.0]) + plt.plot(x1) + plt.plot(x2) + plt.plot(x3) + plt.xlim(0, fs * 0.01) + plt.savefig("q06.pdf") diff --git a/dsugawara/chapter05/q07.py b/dsugawara/chapter05/q07.py new file mode 100644 index 0000000..d2079da --- /dev/null +++ b/dsugawara/chapter05/q07.py @@ -0,0 +1,152 @@ +import numpy as np +import matplotlib.pyplot as plt +import q06 + + +def padding(L, S, x): + """N点の信号xにゼロ詰めを行う + Args: + L(int): 窓幅 + S(int): シフト幅 + x(ndarray): 入力信号 + Return + x_pad: ゼロ詰め後の信号 + """ + + x_pad = np.pad(x, (L - S, L - S)) + temp = S - len(x_pad) % S + x_pad = np.pad(x_pad, (0, temp)) + return x_pad + + +def flame_div(L, S, x): + """N点の信号xをフレーム分割 + Args: + L(int): 窓幅 + S(int): シフト幅 + x(ndarray): 入力信号 + Return + x_div: フレーム分割したx + """ + x_pad = padding(L, S, x) + N = len(x_pad) + T = ((N - L) // S) + 1 + x_div = np.empty([T, L]) + for t in range(0, T): + x_div[t] = x_pad[t * S : t * S + L] + + return x_div + + +def my_stft(L, S, x, win): + """短時間フーリエ変換 + Args: + L(int): 窓幅 + S(int): シフト幅 + x(ndarray): 入力信号 + w(ndarray): 窓関数 + Return + X: 変換後の信号 + """ + x_div = flame_div(L, S, x) + T = x_div.shape[0] + X = np.empty([L // 2 + 1, T], dtype="complex") + for t in range(0, T): + X[:, t] = np.fft.rfft(x_div[t, :] * win) + + return X + + +def optimal_win(S, win): + """合成窓を作成 + Args: + S(int): シフト幅 + win(ndarray): 窓 + Return: + win_opt: 最適合成窓 + """ + + L = win.size + Q = L // S + win_opt = np.empty(L) + for l in range(0, L): + sigma = 0 + for m in range(-(Q - 1), Q): + if l - m * S >= 0 and l - m * S < L: + sigma += win[l - m * S] ** 2 + win_opt[l] = win[l] / sigma + + return win_opt + + +def my_istft(S, X, win): + """逆短時間フーリエ変換 + Args: + S(int): シフト幅 + X(ndarray): 入力信号 + win(ndarray): 窓関数 + Return: + x: 変換後の信号 + """ + F = X.shape[0] + T = X.shape[1] + N = 2 * (F - 1) + M = S * (T - 1) + N + + x = np.zeros(M) + z = np.empty([T, N]) + for t in range(0, T): + z[t, :] = np.fft.irfft(X[:, t]) + x[t * S : t * S + N] = x[t * S : t * S + N] + win * z[t, :] + + return x + + +def delayed_sum_beamformer(L, S, win, X): + _, F, T = X.shape + tau = np.array([0, 10 / fs, 20 / fs]) + w = np.empty([F, 3], dtype="complex") + Y = np.empty([F, T], dtype="complex") + for Fn in range(0, F): + f = (fs / 2) / (F - 1) * Fn + w[Fn, :] = 1 / 3 * np.exp(-2j * np.pi * f * tau) + Y[Fn, :] = np.dot(np.conjugate(w[Fn, :].T), X[:, Fn, :]) + + y = my_istft(S, Y, win) + + return y + + +if __name__ == "__main__": + # 信号の生成 + A = 1.0 + fs = 16000 + sec = 1 + f = 440 + t = np.arange(0, fs * sec) / fs + + s = A * np.sin(2 * np.pi * f * t) + noise = q06.make_noise(s, 10) + + x1 = s + noise + s2 = np.pad(s, [10, 0]) + x2 = s2[0 : len(s2) - 10] + noise + s3 = np.pad(s, [20, 0]) + x3 = s3[0 : len(s3) - 20] + noise + + L = 1024 + S = 512 + win = np.hanning(L) + X1 = my_stft(L, S, x1, win) + X2 = my_stft(L, S, x2, win) + X3 = my_stft(L, S, x3, win) + X = np.stack([X1, X2, X3]) + + y = delayed_sum_beamformer(L, S, win, X) + x1_pad = padding(L, S, x1) + + plt.figure(figsize=[6.0, 4.0]) + plt.plot(np.arange(0, sec * y.shape[0]) / fs, y) + plt.plot(np.arange(0, sec * y.shape[0]) / fs, x1_pad) + plt.xlim(0.02, 0.04) + plt.savefig("q07.pdf") diff --git a/dsugawara/chapter05/q08.py b/dsugawara/chapter05/q08.py new file mode 100644 index 0000000..19c7041 --- /dev/null +++ b/dsugawara/chapter05/q08.py @@ -0,0 +1,48 @@ +import numpy as np +import matplotlib.pyplot as plt +import q03 + + +def plot_beampattern(w, p, fs, name="q08"): + """ビームパターンを描画 + Args: + w (ndarray): 周波数領域におけるビームフォーマのフィルタ + p (ndarray): マイクアレイの座標 + fs (int): サンプリング周波数 + """ + F = w.shape[0] + M = p.shape[1] + a = np.empty([M, F, 360], dtype="complex") + psi = np.empty([F, 360], dtype="complex") + for theta in range(0, 360): + for f in range(0, F): + a[:, f, theta] = q03.array_vector(p, theta, (fs / 2) / (F - 1) * f) + psi[f, theta] = np.dot(np.conjugate(w[f, :].T), a[:, f, theta]) + + plt.figure(figsize=[6.0, 4.0]) + plt.pcolormesh( + np.arange(0, 360), + np.arange(F), + 20 * np.log10(np.abs(psi)), + cmap="gray", + vmin=-20, + ) + plt.colorbar() + plt.savefig(f"{name}_d={p[2,0]}.pdf") + + +if __name__ == "__main__": + L = 1024 + S = 512 + F = L // 2 + 1 + fs = 16000 + d = 0.05 + p = np.array([[-d, 0, 0], [0, 0, 0], [d, 0, 0]]) + + w = np.empty([F, 3], dtype="complex") + for Fn in range(0, F): + f = (fs / 2) / (F - 1) * Fn + w[Fn, :] = q03.array_vector(p, 0, f) + w = w / 3 + + plot_beampattern(w, p, fs) diff --git a/dsugawara/chapter05/q09.py b/dsugawara/chapter05/q09.py new file mode 100644 index 0000000..9961da4 --- /dev/null +++ b/dsugawara/chapter05/q09.py @@ -0,0 +1,21 @@ +import numpy as np +import matplotlib.pyplot as plt +import q03 +import q08 + +if __name__ == "__main__": + L = 1024 + S = 512 + F = L // 2 + 1 + fs = 16000 + + d_list = [2, 5, 10] # マイク間隔 + for d in d_list: + p = np.array([[-d, 0, 0], [0, 0, 0], [d, 0, 0]]) + w = np.empty([F, 3], dtype="complex") + for Fn in range(0, F): + f = (fs / 2) / (F - 1) * Fn + w[Fn, :] = q03.array_vector(p, 0, f) + w = w / 3 + + q08.plot_beampattern(w, p, fs, name="q09") diff --git a/dsugawara/chapter05/q10.py b/dsugawara/chapter05/q10.py new file mode 100644 index 0000000..3609ccf --- /dev/null +++ b/dsugawara/chapter05/q10.py @@ -0,0 +1,67 @@ +import numpy as np +import matplotlib.pyplot as plt +import q01 +import q03 +import q05 +import q06 +import q07 +import q08 + + +def spatial_spectrum(z): + + M = z.shape[0] + + L = 1024 + S = 512 + F = L // 2 + 1 + win = np.hanning(L) + + Z = [] + for m in range(0, M): + Zm = q07.my_stft(L, S, z[m, :], win) + Z.append(Zm) + Z = np.array(Z) + + R = q05.scm(Z) + + F = Z.shape[1] + w = np.empty([F, 360, M], dtype="complex") + for f in range(0, F): + for theta in range(0, 360): + w[f, theta, :] = q01.linear_array_vector(0.05, M, theta, f) + + P = np.empty([F, 360], dtype="complex") + for f in range(0, F): + for theta in range(0, 360): + P[f, theta] = np.conjugate(w[f, theta, :]).T @ R[f, :, :] @ w[f, theta, :] + + for f in range(21, 30): + plt.figure(figsize=[6.0, 4.0]) + plt.plot(np.arange(0, 360), 20 * np.log10(np.abs(P[f, :]))) + plt.savefig(f"q10_f={f}.pdf") + + +if __name__ == "__main__": + M = 3 + d = 0.05 + + A = 1.0 + fs = 16000 + sec = 1 + f = 440 + t = np.arange(0, fs * sec) / fs + + s = A * np.sin(2 * np.pi * f * t) + + noise = q06.make_noise(s, 10) + + x1 = s + noise + s2 = np.pad(s, [10, 0]) + x2 = s2[0 : len(s2) - 10] + noise + s3 = np.pad(s, [20, 0]) + x3 = s3[0 : len(s3) - 20] + noise + + z = np.stack([x1, x2, x3]) + + spatial_spectrum(z)