Merge pull request #22 from kuu8902/kkoiso/chapter04

４章を追加 koiso

Author: taishi-n

kkoiso/chapter04/q01.py

@@ -0,0 +1,12 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def zero_padding(x, L, S):
+    N = len(x)
+    padding_size = L - S
+    padded_x = np.concatenate(([0] * padding_size, x, [0] * padding_size))
+    zero_add = S - (len(padded_x) % S) + (L - S)
+    padded_x = np.concatenate((padded_x, [0] * zero_add))
+
+    return padded_x

kkoiso/chapter04/q02.py

@@ -0,0 +1,14 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q01 import zero_padding
+
+
+def frame_split(x, L, S):
+    padded_x = zero_padding(x, L, S)
+    T = (len(padded_x) - L) // S + 1
+
+    frames = np.zeros((T, L))
+    for t in range(T):
+        frames[t] = padded_x[t * S : t * S + L]
+
+    return frames

kkoiso/chapter04/q03.py

@@ -0,0 +1,11 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q02 import frame_split
+
+
+def stft(x, L, S, w):
+    frames = frame_split(x, L, S)
+    windowed_frames = frames * w
+    stft_matrix = np.fft.rfft(windowed_frames, axis=1)
+
+    return stft_matrix

kkoiso/chapter04/q04.py

@@ -0,0 +1,40 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import stft
+
+
+fs = 16000
+f = 440
+a = 1.0
+t = np.arange(0, 0.1, 1 / fs)
+x = a * np.sin(2 * np.pi * f * t)
+
+
+L = 1000
+S = 500
+w = np.hamming(L)
+
+stft_matrix = stft(x, L, S, w)
+
+Amplitude_spectrogram = np.abs(stft_matrix)
+phase_spectrogram = np.angle(stft_matrix)
+
+num_frames, num_freq_bins = Amplitude_spectrogram.shape
+time_axis = np.arange(num_frames) * (S / fs)
+freq_axis = np.linspace(0, fs / 2, num_freq_bins)
+plt.figure(figsize=(10, 8))
+
+plt.subplot(2, 1, 1)
+plt.title("Ampliitude Spectrogram")
+plt.pcolormesh(time_axis, freq_axis, Amplitude_spectrogram.T)
+plt.ylabel("Frequency")
+plt.xlabel("Frame")
+
+plt.subplot(2, 1, 2)
+plt.title("Phase Spectrogram")
+plt.pcolormesh(time_axis, freq_axis, phase_spectrogram.T)
+plt.ylabel("Frequency")
+plt.xlabel("Frame")
+
+plt.tight_layout()
+plt.show()

kkoiso/chapter04/q05.py

@@ -0,0 +1,28 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+import numpy as np
+
+
+def synthesis_window(w, S):
+    L = len(w)
+    Q = L // S
+    w_s = np.zeros(L)
+
+    for l in range(L):
+        sum_w = 0
+        for m in range(-(Q - 1), Q):
+            idx = l - m * S
+            if 0 <= idx < L:
+                sum_w += w[idx] ** 2
+        if sum_w > 0:
+            w_s[l] = w[l] / sum_w
+        else:
+            w_s[l] = 0
+
+    return w_s
+
+
+w = np.ones(6)
+print(synthesis_window(w, 2))

kkoiso/chapter04/q06.py

@@ -0,0 +1,22 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q05 import synthesis_window
+
+import numpy as np
+
+
+def istft(X, S, w):
+
+    F, T = X.shape
+    N = 2 * (F - 1)
+    M = S * (T - 1) + N
+
+    x_hat = np.zeros(M)
+    w_s = synthesis_window(w, S)
+    for t in range(T):
+        z_t = np.fft.irfft(X[:, t], n=N)
+
+        for n in range(N):
+            x_hat[t * S + n] += w_s[n] * z_t[n]
+
+    return x_hat

kkoiso/chapter04/q07.py

@@ -0,0 +1,31 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import stft
+from q06 import istft
+from q05 import synthesis_window
+
+fs = 16000
+f = 440
+a = 1.0
+t = np.arange(0, 0.1, 1 / fs)
+x = a * np.sin(2 * np.pi * f * t)
+
+
+L = 1000
+S = 500
+w = np.hamming(L)
+w_s = synthesis_window(w, S)
+stft_matrix = stft(x, L, S, w).T
+x_reconstructed = istft(stft_matrix, S, w)
+plt.figure(figsize=(10, 4))
+
+plt.subplot(2, 1, 1)
+plt.plot(x)
+plt.title("Original Signal")
+
+plt.subplot(2, 1, 2)
+plt.plot(x_reconstructed)
+plt.title("Reconstructed Signal")
+plt.grid(True)
+plt.tight_layout()
+plt.show()

kkoiso/chapter04/q08.py

@@ -0,0 +1,48 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import stft
+
+
+def unity_synthesis_window(X, S):
+    F, T = X.shape
+    N = 2 * (F - 1)
+    M = S * (T - 1) + N
+
+    x_hat = np.zeros(M)
+
+    for t in range(T):
+        z_t = np.fft.irfft(X[:, t], n=N)
+
+        for n in range(N):
+            x_hat[t * S + n] += z_t[n]
+
+    return x_hat
+
+
+fs = 16000
+f = 440
+a = 1.0
+t = np.arange(0, 0.1, 1 / fs)
+x = a * np.sin(2 * np.pi * f * t)
+
+
+L = 1000
+S = 500
+w = np.hamming(L)
+stft_matrix = stft(x, L, S, w).T
+
+x_reconstructed_unity_ws = unity_synthesis_window(stft_matrix, S)
+
+
+plt.figure(figsize=(10, 4))
+
+plt.subplot(2, 1, 1)
+plt.plot(x)
+plt.title("Original Signal")
+
+plt.subplot(2, 1, 2)
+plt.plot(x_reconstructed_unity_ws)
+plt.title("Reconstructed Signal with Unity Synthesis Window")
+plt.grid(True)
+plt.tight_layout()
+plt.show()

kkoiso/chapter04/q09.py

@@ -0,0 +1,27 @@
+from scipy.signal import chirp
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import stft
+
+fs = 16000
+t = np.linspace(0, 1, fs)
+chirp_signal = chirp(t, f0=100, f1=16000, t1=1, method="linear")
+
+
+params = [(100, 50), (200, 100), (400, 200), (800, 400)]
+
+
+plt.figure(figsize=(10, 8))
+for i, (L, S) in enumerate(params):
+    w = np.hamming(L)
+    stft_matrix = stft(chirp_signal, L, S, w)
+    Amplitude_spectrogram = 2 * np.log(np.abs(stft_matrix))
+
+    plt.subplot(len(params), 1, i + 1)
+    plt.pcolormesh(Amplitude_spectrogram.T)
+    plt.title(f"Spectrogram Window Size {L} and Shift {S}")
+    plt.ylabel("Frequency Bin")
+    plt.xlabel("Frame")
+
+plt.tight_layout()
+plt.show()

kkoiso/chapter04/q10.py

@@ -0,0 +1,46 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q03 import stft
+
+
+def convert(spectrogram, fs, S):
+    num_frames, num_freq_bins = spectrogram.shape
+    time_axis = np.arange(num_frames) * (S / fs)
+    freq_axis = np.linspace(0, fs / 2, num_freq_bins)
+    return time_axis, freq_axis
+
+
+fs = 16000
+f = 440
+a = 1.0
+t = np.arange(0, 0.1, 1 / fs)
+x = a * np.sin(2 * np.pi * f * t)
+
+
+L = 1000
+S = 500
+w = np.hamming(L)
+
+stft_matrix = stft(x, L, S, w)
+
+Amplitude_spectrogram = np.abs(stft_matrix)
+phase_spectrogram = np.angle(stft_matrix)
+
+time_axis, freq_axis = convert(Amplitude_spectrogram, fs, S)
+
+plt.figure(figsize=(10, 8))
+
+plt.subplot(2, 1, 1)
+plt.title("Ampliitude Spectrogram")
+plt.pcolormesh(time_axis, freq_axis, Amplitude_spectrogram.T)
+plt.ylabel("Frequency")
+plt.xlabel("Frame")
+
+plt.subplot(2, 1, 2)
+plt.title("Phase Spectrogram")
+plt.pcolormesh(time_axis, freq_axis, phase_spectrogram.T)
+plt.ylabel("Frequency")
+plt.xlabel("Frame")
+
+plt.tight_layout()
+plt.show()