Merge pull request #25 from woodnx/skotsugi/chapter04

Add chapter04

Author: taishi-n

skotsugi/chapter04/q01.py

@@ -0,0 +1,24 @@
+import numpy as np
+
+def zero_padding(L: int, S: int, x):
+  # エラー処理
+  if L < S:
+    raise Exception('Shift size must be smaller than window size.')
+
+  # 出力信号長の算出
+  N = len(x)
+  M = L - S    # > 0
+  A = N + M
+  B = A // S + 1 # 条件を満たす最小のSの倍数
+  C = S * B + M
+
+  # 出力信号
+  tilde = np.zeros(C) # 全体をゼロ埋め
+  tilde[M:M+N] = x
+
+  return tilde
+
+if __name__ == "__main__":
+  x = np.ones(5)
+  print(x)
+  print(zero_padding(4, 3,))

skotsugi/chapter04/q02.py

@@ -0,0 +1,19 @@
+import numpy as np
+from q01 import zero_padding
+
+def divide_frame(L: int, S: int, x):
+  # エラー処理
+  if L < S:
+    raise Exception('Shift size must be smaller than window size.')
+  
+  x_tilde = zero_padding(L, S, x)
+  N = len(x_tilde)
+
+  T = np.arange((N - L) // S + 1)
+  return np.array([ x_tilde[t*S:t*S+L] for t in T ])
+
+if __name__ == "__main__":
+  x = np.ones(5)
+  print(x)
+  print(zero_padding(4, 3, x))
+  print(divide_frame(4, 3, x))

skotsugi/chapter04/q03.py

@@ -0,0 +1,13 @@
+import numpy as np
+from q02 import divide_frame
+
+def stft(L: int, S: int, w, x):
+  # エラー処理
+  if L < S:
+    raise Exception('Shift size must be smaller than window size.')
+
+  flames = divide_frame(L, S, x)
+
+  ffted = np.fft.rfft(flames*w).T
+
+  return ffted

skotsugi/chapter04/q04.png

@@ -0,0 +1,47 @@
+import numpy as np
+
+def sin(A, f, fs, sec):
+  n = np.arange(0, fs*sec) / fs
+  return A * np.sin(2 * np.pi * f * n)
+
+def hamming(N: int): 
+  n = np.arange(N)
+  return 0.54 - 0.46 * np.cos(2 * np.pi * n / (N - 1))
+
+if __name__ == "__main__":
+  import matplotlib.pyplot as plt
+  from q03 import stft
+  
+  L = 1000
+  S = 500
+  sec = 0.1
+  freq = 440
+  fs = 16000
+  x = sin(1, freq, fs, sec)
+  w = hamming(L)
+
+  X = stft(L, S, w, x)
+  F = X.shape[0]
+  T = X.shape[1]
+
+  t = (np.arange(T) + 0.5) * sec / T
+  f = np.arange(F) / F * fs / 2
+
+  A = 20*np.log10(np.abs(X))
+  ang = np.angle(X)
+
+  fig, ax = plt.subplots(1, 2)
+
+  cb = ax[0].pcolormesh(t, f, A)
+  ax[0].set_ylabel('Frequency [Hz]')
+  ax[0].set_xlabel('Time [sec]')
+
+  cb = ax[1].pcolormesh(t, f, ang)
+  ax[1].set_ylabel('Frequency [Hz]')
+  ax[1].set_xlabel('Time [sec]')
+
+  fig.colorbar(cb, orientation="vertical", ax=ax)
+
+  plt.grid()
+  plt.tight_layout()
+  plt.savefig('skotsugi/chapter04/q04.png')

skotsugi/chapter04/q04.py

@@ -0,0 +1,25 @@
+import numpy as np
+
+def synthesis_window(S: int, w: np.ndarray): 
+  L = w.size
+  Q = L // S
+  
+  s = np.array([
+    np.sum([
+      w[l - m*S]**2 
+      for m in range(-(Q-1), Q) 
+      if 0 <= l - m*S and L > l - m*S
+    ]) for l in range(L) 
+  ])
+
+  return w / s
+
+
+if __name__ == "__main__":
+  import matplotlib.pyplot as plt
+  
+  N = 500
+  hamming = np.hamming(N)
+  ws = synthesis_window(N, hamming)
+  plt.plot(ws)
+  plt.savefig("./skotsugi/chapter04/q05.png")

skotsugi/chapter04/q05.py

@@ -0,0 +1,24 @@
+import numpy as np
+from q05 import synthesis_window
+
+def istft(S: int, w: np.ndarray, X: np.ndarray):
+  F, T = X.shape
+  N = 2 * (F - 1)
+  M = S * (T - 1) + N
+  
+  hat = np.zeros(M)
+  
+  z = np.array([ np.fft.irfft(X[:,t]) for t in range(T) ])
+
+  if z.shape[1] != N:
+    raise "len(z[n]) must be 2(F - 1)"
+  
+  ws = synthesis_window(S, w)
+
+  for t in range(T):
+    idx = t * S
+
+    # Overlap add
+    hat[idx:idx + N] += ws * z[t]
+
+  return hat

skotsugi/chapter04/q06.py

@@ -0,0 +1,20 @@
+import matplotlib.pyplot as plt
+from q03 import stft
+from q04 import sin, hamming
+from q06 import istft
+
+L = 1000
+S = 500
+sec = 0.1
+freq = 440
+fs = 16000
+x = sin(1, freq, fs, sec)
+w = hamming(L)
+
+X = stft(L, S, w, x)
+x_istft = istft(S, w, X)
+
+plt.plot(x_istft)
+plt.grid()
+plt.xlim(0)
+plt.savefig('skotsugi/chapter04/q07.png')

skotsugi/chapter04/q07.png

@@ -0,0 +1,43 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from q03 import stft
+from q04 import sin, hamming
+
+def istft(S: int, X: np.ndarray):
+  F, T = X.shape
+  N = 2 * (F - 1)
+  M = S * (T - 1) + N
+  
+  hat = np.zeros(M)
+  
+  z = np.array([ np.fft.irfft(X[0:F,t]) for t in range(T) ])
+
+  if z.shape[1] != N:
+    raise "len(z[n]) must be 2(F - 1)"
+  
+  ws = np.ones(N)
+
+  for t in range(T):
+    idx = t * S
+
+    # Overlap add
+    hat[idx:idx + N] = hat[idx:idx + N] + ws * z[t]
+
+  return hat
+
+
+L = 1000
+S = 500
+sec = 0.1
+freq = 440
+fs = 16000
+x = sin(1, freq, fs, sec)
+w = hamming(L)
+
+X = stft(L, S, w, x)
+x_istft = istft(S, X)
+
+plt.plot(x_istft)
+plt.grid()
+plt.xlim(0)
+plt.savefig('skotsugi/chapter04/q08.png')

skotsugi/chapter04/q07.py

@@ -0,0 +1,35 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.signal import chirp
+from q03 import stft
+from q04 import hamming
+
+sec = 1
+fs = 16000
+t = np.linspace(0, sec, fs+1) 
+x = chirp(t, f0=100, f1=16000, t1=sec)
+
+L = 50
+S = 25
+Xs = []
+for i in range(1, 5):
+  L = L * 2
+  S = S * 2
+  w = hamming(L)
+  X = stft(L, S, w, x)
+  Xs.append(X)
+
+fig, ax = plt.subplots(4, 1, sharey=True, sharex=True)
+for i, X in enumerate(Xs):
+  A = np.abs(X)
+
+  F = X.shape[0]
+  T = X.shape[1]
+
+  cb = ax[i].pcolormesh(A)
+
+plt.tight_layout()
+
+fig.colorbar(cb, ax=ax)
+
+plt.savefig('skotsugi/chapter04/q09.png')

skotsugi/chapter04/q08.png

@@ -0,0 +1,50 @@
+import numpy as np
+
+def to_quantity(X: np.ndarray, fs: int, S: int):
+  F = X.shape[0] # 1フレームあたりの信号長
+  T = X.shape[1] # フレーム数
+  L = 2 * (F - 1)
+
+  t = np.arange(T) * S / fs
+  f = np.arange(F) * fs / L
+
+  print(t)
+
+  return t, f
+
+if __name__ == "__main__":
+  import matplotlib.pyplot as plt
+  from q03 import stft
+  from q04 import sin, hamming
+  
+  L = 1000
+  S = 500
+  sec = 0.1
+  freq = 440
+  fs = 16000
+  x = sin(1, freq, fs, sec)
+  w = hamming(L)
+
+  X = stft(L, S, w, x)
+
+  t, f = to_quantity(X, fs, S)
+
+  A = 20 * np.log10(np.abs(X))
+  ang = np.angle(X)
+
+  fig, ax = plt.subplots(1, 2)
+
+  cb = ax[0].pcolormesh(t, f, A)
+  ax[0].set_ylabel('Frequency [Hz]')
+  ax[0].set_xlabel('Time [sec]')
+
+  cb = ax[1].pcolormesh(t, f, ang)
+  ax[1].set_ylabel('Frequency [Hz]')
+  ax[1].set_xlabel('Time [sec]')
+
+  plt.tight_layout()
+
+  fig.colorbar(cb, ax=ax)
+
+  plt.grid()
+  plt.savefig('skotsugi/chapter04/q10.png')