add chapter04

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)

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)

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)

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()

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()

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

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()

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()