diff --git a/kkazama/chapter04/q01.py b/kkazama/chapter04/q01.py new file mode 100644 index 0000000..2a903f3 --- /dev/null +++ b/kkazama/chapter04/q01.py @@ -0,0 +1,20 @@ +import numpy as np + + +# ゼロ詰め関数 +def zero_pad(x, L, S): + x = np.pad(x, L - S) + + if len(x) % S != 0: + x = np.pad(x, (0, S - len(x) % S)) + + return x + + +# 確認 +x = [1, 2, 3] +L = 10 +S = 5 +x_pad = zero_pad(x, L, S) + +print(x_pad) diff --git a/kkazama/chapter04/q02.py b/kkazama/chapter04/q02.py new file mode 100644 index 0000000..04a1704 --- /dev/null +++ b/kkazama/chapter04/q02.py @@ -0,0 +1,34 @@ +import numpy as np + + +# ゼロ詰め関数 +def zero_pad(x, L, S): + x = np.pad(x, L - S) + + if len(x) % S != 0: + x = np.pad(x, (0, S - len(x) % S)) + + return x + + +# フレーム分割 +def frame_division(x, L, S): + x = zero_pad(x, L, S) + T = (len(x) - L) // S + 1 + + x_frame = np.zeros((T, L)) + + for t in range(T): + x_frame[t] = x[t * S : t * S + L] + + return x_frame.T + + +# 確認 +x = [1, 1, 1, 1, 1, 1, 1, 1] +L = 4 +S = 3 +x_frame = frame_division(x, L, S) + +print(x_frame) +print(x_frame.shape) diff --git a/kkazama/chapter04/q03.py b/kkazama/chapter04/q03.py new file mode 100644 index 0000000..b7a0f39 --- /dev/null +++ b/kkazama/chapter04/q03.py @@ -0,0 +1,47 @@ +import numpy as np + + +# ゼロ詰め関数 +def zero_pad(x, L, S): + x = np.pad(x, L - S) + + if len(x) % S != 0: + x = np.pad(x, (0, S - len(x) % S)) + + return x + + +# フレーム分割 +def frame_division(x, L, S): + x = zero_pad(x, L, S) + T = (len(x) - L) // S + 1 + + x_frame = np.zeros((T, L)) + + for t in range(T): + x_frame[t] = x[t * S : t * S + L] + + return x_frame.T + + +# 短時間フーリエ変換 +def stft(x, L, S, w): + x_frame = frame_division(x, L, S).T + T = len(x_frame) + X = np.zeros((T, L // 2 + 1), dtype=complex) + + for t in range(T): + x_t = x_frame[t] * w + X[t] = np.fft.rfft(x_t) + + return X.T + + +# 確認 +L = 4 +S = 3 +x = np.ones(8) +wnd = np.hamming(L) +x_stft = stft(x, L, S, wnd) + +print(x_stft) diff --git a/kkazama/chapter04/q04.py b/kkazama/chapter04/q04.py new file mode 100644 index 0000000..a6332e2 --- /dev/null +++ b/kkazama/chapter04/q04.py @@ -0,0 +1,67 @@ +import numpy as np +import matplotlib.pyplot as plt + + +# ゼロ詰め関数 +def zero_pad(x, L, S): + x = np.pad(x, L - S) + + if len(x) % S != 0: + x = np.pad(x, (0, S - len(x) % S)) + + return x + + +# フレーム分割 +def frame_division(x, L, S): + x = zero_pad(x, L, S) + T = (len(x) - L) // S + 1 + + x_frame = np.zeros((T, L)) + + for t in range(T): + x_frame[t] = x[t * S : t * S + L] + + return x_frame.T + + +# 短時間フーリエ変換 +def stft(x, L, S, w): + x_frame = frame_division(x, L, S).T + T = len(x_frame) + X = np.zeros((T, L // 2 + 1), dtype=complex) + + for t in range(T): + x_t = x_frame[t] * w + X[t] = np.fft.rfft(x_t) + + return X.T + + +# 確認 +A = 1 +f = 440 +fs = 16000 +sec = 0.1 + +t = np.arange(fs * sec) / fs +x = A * np.sin(2 * np.pi * f * t) + +L = 1000 +S = 500 +w = np.hamming(L) + +X = stft(x, L, S, w) + +# プロット +fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(6, 5)) + +p0 = ax[0].pcolormesh(np.abs(X)) +ax[0].set_title("Amplitude spectrum") + +p1 = ax[1].pcolormesh(np.angle(X)) +ax[1].set_title("Phase spectrum") + +plt.tight_layout() + +plt.show() diff --git a/kkazama/chapter04/q05.py b/kkazama/chapter04/q05.py new file mode 100644 index 0000000..bb434ef --- /dev/null +++ b/kkazama/chapter04/q05.py @@ -0,0 +1,35 @@ +import numpy as np +import matplotlib.pyplot as plt + + +# 合成窓 +def window_synth(w, L, S): + Q = L // S + w_synth = np.zeros(L) + + for l in range(L): + w_shift_sum = 0 + for m in range(-(Q - 1), Q): + if (l - m * S) >= 0 and (l - m * S) < L: + w_shift_sum += w[l - m * S]**2 + + w_synth[l] = w[l] / w_shift_sum + + return w_synth + + +# 確認 +L = 1000 +S = 250 +w = np.hamming(L) + +w_synth = window_synth(w, L, S) + +a = np.sum(w * w_synth) +print(a) + +plt.subplot(2, 1, 1) +plt.plot(w) +plt.subplot(2, 1, 2) +plt.plot(w_synth) +plt.show() diff --git a/kkazama/chapter04/q06.py b/kkazama/chapter04/q06.py new file mode 100644 index 0000000..9fa8a36 --- /dev/null +++ b/kkazama/chapter04/q06.py @@ -0,0 +1,38 @@ +import numpy as np +import matplotlib.pyplot as plt + + +# 合成窓 +def window_synth(w, L, S): + Q = L // S + w_synth = np.zeros(L) + + for l in range(L): + w_shift_sum = 0 + for m in range(-(Q - 1), Q): + if (l - m * S) >= 0 and (l - m * S) < L: + w_shift_sum += w[l - m * S]**2 + + w_synth[l] = w[l] / w_shift_sum + + return w_synth + + +def istft(X, S): + F = X.shape[0] + T = X.shape[1] + + N = 2 * (F - 1) + M = S * (T - 1) + N + + x = np.zeros(M) + z = np.zeros((T, N)) + w = np.hamming(N) + w_synth = window_synth(w) + + for t in range(T): + for n in range(N): + z[t][n] = np.fft.irfft(X[:, t])[n] + x[t * S + n] = x[t * S + n] + w_synth[n] * z[t][n] + + return x diff --git a/kkazama/chapter04/q07.py b/kkazama/chapter04/q07.py new file mode 100644 index 0000000..d4180d2 --- /dev/null +++ b/kkazama/chapter04/q07.py @@ -0,0 +1,61 @@ +import numpy as np +import matplotlib.pyplot as plt +from q04 import stft + + +# 合成窓 +def window_synth(w, L, S): + Q = L // S + w_synth = np.zeros(L) + + for l in range(L): + w_shift_sum = 0 + for m in range(-(Q - 1), Q): + if (l - m * S) >= 0 and (l - m * S) < L: + w_shift_sum += w[l - m * S]**2 + + w_synth[l] = w[l] / w_shift_sum + + return w_synth + + +def istft(X, S): + F = X.shape[0] + T = X.shape[1] + + N = 2 * (F - 1) + M = S * (T - 1) + N + + x = np.zeros(M) + z = np.zeros((T, N)) + w = np.hamming(N) + w_synth = window_synth(w, N, S) + + for t in range(T): + for n in range(N): + z[t][n] = np.fft.irfft(X[:, t])[n] + x[t * S + n] = x[t * S + n] + w_synth[n] * z[t][n] + + return x + + +# 確認 +A = 1 +f = 440 +fs = 16000 +sec = 0.1 + +t = np.arange(fs * sec) / fs +x = A * np.sin(2 * np.pi * f * t) + +L = 1000 +S = 500 +w = np.hamming(L) + +X = stft(x, L, S, w) + +x_istft = istft(X, S) + +plt.plot(x_istft) +plt.xlim([500, 1000]) +plt.show() diff --git a/kkazama/chapter04/q08.py b/kkazama/chapter04/q08.py new file mode 100644 index 0000000..02e1705 --- /dev/null +++ b/kkazama/chapter04/q08.py @@ -0,0 +1,76 @@ +import numpy as np +import matplotlib.pyplot as plt +from q04 import stft +from q07 import istft +from q01 import zero_pad + + +# 合成窓 +def window_synth(w, L, S): + Q = L // S + w_synth = np.zeros(L) + + for l in range(L): + w_shift_sum = 0 + for m in range(-(Q - 1), Q): + if (l - m * S) >= 0 and (l - m * S) < L: + w_shift_sum += w[l - m * S]**2 + + w_synth[l] = w[l] / w_shift_sum + + return w_synth + + +def istft_w(X, S): + F = X.shape[0] + T = X.shape[1] + + N = 2 * (F - 1) + M = S * (T - 1) + N + + x = np.zeros(M) + z = np.zeros((T, N)) + w = np.zeros(N) + 1.0 + + for t in range(T): + for n in range(N): + z[t][n] = np.fft.irfft(X[:, t])[n] + x[t * S + n] = x[t * S + n] + w[n] * z[t][n] + + return x + + +# 確認 +A = 1 +f = 440 +fs = 16000 +sec = 0.1 + +t = np.arange(fs * sec) / fs +x = A * np.sin(2 * np.pi * f * t) + +L = 1000 +S = 500 +w = np.hamming(L) + +X = stft(x, L, S, w) + +x_istft = istft(X, S) +x_istft_w = istft_w(X, S) + +x_pad = zero_pad(x, L, S) + +# 再構成誤差 +diff = np.sum(x_pad - x_istft) ** 2 +diff_w = np.sum(x_pad - x_istft_w) ** 2 + +print(diff) +print(diff_w) + +plt.subplot(2, 1, 1) +plt.plot(x_istft) +plt.xlim([500, 1000]) +plt.subplot(2, 1, 2) +plt.plot(x_istft_w) +plt.xlim([500, 1000]) +plt.show() diff --git a/kkazama/chapter04/q09.py b/kkazama/chapter04/q09.py new file mode 100644 index 0000000..422c9b2 --- /dev/null +++ b/kkazama/chapter04/q09.py @@ -0,0 +1,32 @@ +from scipy.signal import chirp +import numpy as np +import matplotlib.pyplot as plt +from q04 import stft + +# チャープ信号 +fs = 16000 +sec = 1 +t = np.linspace(0, 1, fs) +y = chirp(t, f0=100, f1=16000, t1=sec, method='linear') + +L0 = 100 +S0 = 50 +Y_amp = [] + +for i in range(1, 5): + L = i * L0 + S = i * S0 + w = np.hamming(L) + Y = stft(y, L, S, w) + Y_amp.append(np.abs(Y)) + +# プロット +fig, ax = plt.subplots(nrows=4, ncols=1, figsize=(6, 10)) + +for i in range(4): + p = ax[i].pcolormesh(Y_amp[i]) + ax[i].set_title("Amplitude spectrum") + +plt.tight_layout() + +plt.show() diff --git a/kkazama/chapter04/q10.py b/kkazama/chapter04/q10.py new file mode 100644 index 0000000..c96d651 --- /dev/null +++ b/kkazama/chapter04/q10.py @@ -0,0 +1,79 @@ +import numpy as np +import matplotlib.pyplot as plt + + +def calc_axis(spec, fs, S): + F = spec.shape[0] + T = spec.shape[1] + L = (F - 1) * 2 + + F_axis = np.arange(0, F + 1) / F * (fs / 2) + T_axis = (np.arange(0, T + 1) * S - (L - S)) / fs + + return F_axis, T_axis + + +# ゼロ詰め関数 +def zero_pad(x, L, S): + x = np.pad(x, L - S) + + if len(x) % S != 0: + x = np.pad(x, (0, S - len(x) % S)) + + return x + + +# フレーム分割 +def frame_division(x, L, S): + x = zero_pad(x, L, S) + T = (len(x) - L) // S + 1 + + x_frame = np.zeros((T, L)) + + for t in range(T): + x_frame[t] = x[t * S : t * S + L] + + return x_frame.T + + +# 短時間フーリエ変換 +def stft(x, L, S, w): + x_frame = frame_division(x, L, S).T + T = len(x_frame) + X = np.zeros((T, L // 2 + 1), dtype=complex) + + for t in range(T): + x_t = x_frame[t] * w + X[t] = np.fft.rfft(x_t) + + return X.T + + +# 確認 +A = 1 +f = 440 +fs = 16000 +sec = 0.1 + +t = np.arange(fs * sec) / fs +x = A * np.sin(2 * np.pi * f * t) + +L = 1000 +S = 500 +w = np.hamming(L) + +X = stft(x, L, S, w) +F, T = calc_axis(X, fs, S) + +# プロット +fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(6, 5)) + +p0 = ax[0].pcolormesh(T, F, np.abs(X)) +ax[0].set_title("Amplitude spectrum") + +p1 = ax[1].pcolormesh(T, F, np.angle(X)) +ax[1].set_title("Phase spectrum") + +plt.tight_layout() + +plt.show()