diff --git a/eimamura/chapter05/q06.py b/eimamura/chapter05/q06.py new file mode 100644 index 0000000..61276a8 --- /dev/null +++ b/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() diff --git a/eimamura/chapter05/q07.py b/eimamura/chapter05/q07.py new file mode 100644 index 0000000..fcb0e7e --- /dev/null +++ b/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() diff --git a/eimamura/chapter05/q08.py b/eimamura/chapter05/q08.py new file mode 100644 index 0000000..ef49891 --- /dev/null +++ b/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) diff --git a/eimamura/chapter05/q09.py b/eimamura/chapter05/q09.py new file mode 100644 index 0000000..18b5b1b --- /dev/null +++ b/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) diff --git a/eimamura/chapter05/q10.py b/eimamura/chapter05/q10.py new file mode 100644 index 0000000..74ce518 --- /dev/null +++ b/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()