４章後半を追加

Author: 153hashimoto

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,46 @@
+# 時間周波数のインデクスと物理量との対応
+
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import stft
+
+
+def fr_convert(x, fs, S):
+    L = (len(x) - 1) * 2
+    f = np.arange(len(x)) * fs / L
+    return f
+
+
+def t_convert(x, fs, S):
+    t = np.arange(x.shape[1]) * S / fs
+    return 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 = fr_convert(x_stft.T, fs, S)
+t = 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()