Merge pull request #24 from 153hashimoto/khashimoto/chapter04

４章を追加　Hashimoto

Author: taishi-n

khashimoto/chapter04/q01.py

@@ -0,0 +1,31 @@
+# 零詰め
+
+import numpy as np
+
+
+def zeropad(L, S, x):
+    zero1 = np.zeros(L - S)
+    # x_pad = np.concatenate([zero1, x, zero1])  # 1.
+    x_pad = np.concatenate([zero1, x])
+    zero2 = x_pad.size % S
+    if zero2 != 0:
+        x_pad = np.concatenate([x_pad, np.zeros(S - zero2)])  # 2.
+    x_pad = np.concatenate([x_pad, np.zeros(L - S)])
+    return x_pad
+
+
+"""
+# 確認
+x = np.ones(5)
+print(x)
+x_pad = zeropad(4, 2, x)
+print(x_pad)
+
+
+# [1. 1. 1. 1. 1.]
+# [0. 1. 1. 1. 1. 1. 0. 0. 0.]
+
+# [1. 1. 1. 1. 1.]
+# [0. 0. 1. 1. 1. 1. 1. 0. 0. 0.]
+# 不足分が出ないように零詰め
+"""

khashimoto/chapter04/q02.py

@@ -0,0 +1,29 @@
+# フレーム分割
+
+import numpy as np
+from q01 import zeropad
+
+
+def frame(L, S, x):
+    x_pad = zeropad(L, S, x)
+    T = (x_pad.size - L) // S + 1  # フレーム数
+    x_t = np.zeros((T, L))
+    for t in range(T):
+        for l in range(L):
+            x_t[t][l] = x_pad[t * S + l]
+    return x_t
+
+
+"""
+# 確認
+x = np.ones(5)
+print(x)
+print(zeropad(4, 2, x))
+x_2 = frame(4, 2, x)
+print(x_2)
+
+
+# [1. 1. 1. 1. 1.]
+# [[0. 1. 1. 1.]
+# [1. 1. 1. 0.]]
+"""

khashimoto/chapter04/q03.py

@@ -0,0 +1,13 @@
+# STFTの実装
+
+import numpy as np
+from q02 import frame
+
+
+def stft(L, S, w, x):
+    x_t = frame(L, S, x)
+    T = len(x_t)
+    x_stft = np.zeros((T, L // 2 + 1), dtype=complex)
+    for t in range(T):
+        x_stft[t] = np.fft.rfft(x_t[t] * w)
+    return x_stft

khashimoto/chapter04/q04.py

@@ -0,0 +1,33 @@
+# STFTの確認
+
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import stft
+
+A = 1
+f = 440
+fs = 16000
+time = np.arange(fs * 0.1) / fs
+x = A * np.sin(2 * np.pi * f * time)  # 正弦波
+
+L = 1000
+S = 500
+w = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(L) / (L - 1))  # Hamming窓
+x_stft = stft(L, S, w, x)
+
+
+# プロット
+T = np.arange(len(x_stft))  # 時間軸
+fr = np.arange(x_stft.shape[1]) / x_stft.shape[1] * fs / 2  # 周波数軸
+
+plt.subplot(1, 2, 1)
+plt.pcolormesh(T, fr, abs(x_stft).T)
+plt.xlabel("Time[s]")
+plt.ylabel("Frequency[Hz]")
+plt.title("Amplitude spectrum")
+plt.subplot(1, 2, 2)
+plt.pcolormesh(T, fr, np.angle(x_stft).T)
+plt.xlabel("Time[s]")
+plt.ylabel("Frequency[kHz]")
+plt.title("Phase spectrum")
+plt.show()

khashimoto/chapter04/q05.py

@@ -0,0 +1,29 @@
+# 合成窓
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def window(S, w):
+    L = w.size
+    Q = L // S  # 1.
+    w_s = np.zeros(L)
+    for l in range(L):
+        sum = 0
+        for m in range(-(Q - 1), Q):
+            if 0 <= l - m * S < L:
+                sum += w[l - m * S] ** 2  # 2.
+        w_s[l] = w[l] / sum
+    return w_s
+
+
+"""
+# 確認
+L = 1024
+w = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(L) / (L - 1))  # Hamming窓
+plt.plot(window(256, w))
+plt.show()
+
+
+# [0.33333333 0.33333333 0.33333333 0.33333333 0.33333333 0.33333333]
+"""

khashimoto/chapter04/q06.py

@@ -0,0 +1,30 @@
+# ISTFTの実装
+
+import numpy as np
+from q05 import window
+
+
+def istft(S, X):
+    F = len(X)
+    T = X.shape[1]
+
+    # 1.
+    N = 2 * (F - 1)
+    M = S * (T - 1) + N
+
+    # 2.
+    x = np.zeros(M)
+
+    # 3.
+    z = np.zeros((N, T))
+    for t in range(T):
+        z[:, t] = np.fft.irfft(X[:, t])
+
+    # 4.
+    w = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(N) / (N - 1))  # Hamming窓
+    w_s = window(S, w)
+    n = np.arange(N)
+    for t in range(T):
+        x[t * S + n] = x[t * S + n] + w_s[n] * z[n, t]
+
+    return x

khashimoto/chapter04/q07.py

@@ -0,0 +1,22 @@
+# ISTFTの確認
+
+import numpy as np
+import matplotlib.pyplot as plt
+from q06 import istft
+from q03 import stft
+
+A = 1
+f = 440
+fs = 16000
+time = np.arange(fs * 0.1) / fs
+x = A * np.sin(2 * np.pi * f * time)  # 正弦波
+
+L = 1000
+S = 500
+w = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(L) / (L - 1))  # Hamming窓
+x_stft = stft(L, S, w, x)  # 4. の結果
+
+x_istft = istft(S, x_stft.T)
+
+plt.plot(x_istft)
+plt.show()

khashimoto/chapter04/q08.py

@@ -0,0 +1,46 @@
+# 合成窓の確認
+
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import stft
+
+
+def istft_2(S, X):
+    F = len(X)
+    T = X.shape[1]
+
+    # 1.
+    N = 2 * (F - 1)
+    M = S * (T - 1) + N
+
+    # 2.
+    x = np.zeros(M)
+
+    # 3.
+    z = np.zeros((N, T))
+    for t in range(T):
+        z[:, t] = np.fft.irfft(X[:, t])
+
+    # 4.
+    n = np.arange(N)
+    for t in range(T):
+        x[t * S + n] = x[t * S + n] + 1 * z[n, t]  # w_s[n] = 1
+
+    return x
+
+
+A = 1
+f = 440
+fs = 16000
+time = np.arange(fs * 0.1) / fs
+x = A * np.sin(2 * np.pi * f * time)  # 正弦波
+
+L = 1000
+S = 500
+w = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(L) / (L - 1))  # Hamming窓
+x_stft = stft(L, S, w, x)  # 4. の結果
+
+x_istft = istft_2(S, x_stft.T)
+
+plt.plot(x_istft)
+plt.show()

khashimoto/chapter04/q09.py

@@ -0,0 +1,65 @@
+# 不確定性原理の確認
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.signal import chirp
+from q03 import stft
+
+fs = 16000
+t = np.arange(0, fs) / fs
+sig = chirp(t, f0=100, f1=16000, t1=1)
+
+# スペクトログラムのプロット
+
+# 100/50
+L = 100
+S = 50
+w = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(L) / (L - 1))  # Hamming窓
+sig_stft = stft(L, S, w, sig)
+
+plt.subplot(4, 2, 1)
+plt.pcolormesh(abs(sig_stft).T)
+plt.title("Amplitude spectrum (100/50)")
+plt.subplot(4, 2, 2)
+plt.pcolormesh(np.angle(sig_stft).T)
+plt.title("Phase spectrum (100/50)")
+
+# 200/100
+L = 200
+S = 100
+w = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(L) / (L - 1))  # Hamming窓
+sig_stft = stft(L, S, w, sig)
+
+plt.subplot(4, 2, 3)
+plt.pcolormesh(abs(sig_stft).T)
+plt.title("Amplitude spectrum (200/100)")
+plt.subplot(4, 2, 4)
+plt.pcolormesh(np.angle(sig_stft).T)
+plt.title("Phase spectrum (200/100)")
+
+# 400/200
+L = 400
+S = 200
+w = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(L) / (L - 1))  # Hamming窓
+sig_stft = stft(L, S, w, sig)
+
+plt.subplot(4, 2, 5)
+plt.pcolormesh(abs(sig_stft).T)
+plt.title("Amplitude spectrum (400/200)")
+plt.subplot(4, 2, 6)
+plt.pcolormesh(np.angle(sig_stft).T)
+plt.title("Phase spectrum (400/200)")
+
+# 800/400
+L = 800
+S = 400
+w = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(L) / (L - 1))  # Hamming窓
+sig_stft = stft(L, S, w, sig)
+
+plt.subplot(4, 2, 7)
+plt.pcolormesh(abs(sig_stft).T)
+plt.title("Amplitude spectrum (800/400)")
+plt.subplot(4, 2, 8)
+plt.pcolormesh(np.angle(sig_stft).T)
+plt.title("Phase spectrum (800/400)")
+plt.show()

khashimoto/chapter04/q10.py

@@ -0,0 +1,41 @@
+# 時間周波数のインデクスと物理量との対応
+
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import stft
+
+
+def convert(x, fs, S):
+    L = (len(x) - 1) * 2
+    f = np.arange(len(x)) * fs / L
+    t = np.arange(x.shape[1]) * S / fs
+    return f, t
+
+
+# ４. の結果を再度プロット
+A = 1
+f = 440
+fs = 16000
+time = np.arange(fs * 0.1) / fs
+x = A * np.sin(2 * np.pi * f * time)  # 正弦波
+
+L = 1000
+S = 500
+w = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(L) / (L - 1))  # Hamming窓
+x_stft = stft(L, S, w, x)
+
+
+# プロット
+fr, t = convert(x_stft.T, fs, S)
+
+plt.subplot(1, 2, 1)
+plt.pcolormesh(t, fr, abs(x_stft).T)
+plt.xlabel("Time[s]")
+plt.ylabel("Frequency[Hz]")
+plt.title("Amplitude spectrum")
+plt.subplot(1, 2, 2)
+plt.pcolormesh(t, fr, np.angle(x_stft).T)
+plt.xlabel("Time[s]")
+plt.ylabel("Frequency[Hz]")
+plt.title("Phase spectrum")
+plt.show()