Merge pull request #17 from woodnx/skotsugi/chapter03

Add skotsugi/chapter03

Author: taishi-n

skotsugi/chapter03/q01.py

@@ -0,0 +1,17 @@
+import numpy as np
+
+def linear_conv(x, h):
+  N = len(x)
+  M = 2*N - 1
+  z = np.zeros(M, dtype = 'complex_')
+
+  for i in range(M):
+    for k in range(N):
+      j = i - k
+      
+      if j < 0 or j > N - 1:
+        continue
+      else:
+        z[i] += x[k] * h[j]
+
+  return z

skotsugi/chapter03/q02.py

@@ -0,0 +1,11 @@
+import numpy as np
+
+def circular_conv(x, h):
+  N = len(x)
+  z = np.zeros(N, dtype = 'complex_')
+
+  for i in range(N):
+    for k in range(N):
+      z[i] += x[k] * h[(i - k) % N]
+
+  return z

skotsugi/chapter03/q03.py

@@ -0,0 +1,20 @@
+import numpy as np
+
+def circular_zero_conv(x, h):
+  N = len(x)
+  M = N * 2
+  
+  x_zero = np.zeros(M, dtype = 'complex_')
+  x_zero[:N] = x
+
+  h_zero = np.zeros(M, dtype = 'complex_')
+  h_zero[:N] = h
+
+  z = np.zeros(M, dtype = 'complex_')
+
+  for i in range(M):
+    for k in range(M):
+      j = (i - k) % M
+      z[i] += x_zero[k] * h_zero[j]
+
+  return z

skotsugi/chapter03/q04.png

@@ -0,0 +1,23 @@
+import matplotlib.pyplot as plt
+from q01 import linear_conv
+from q02 import circular_conv
+from q03 import circular_zero_conv
+
+x = [4, 3, 2, 1]
+y = [1, 0, -1, 0]
+
+z_l = linear_conv(x, y)
+z_c = circular_conv(x, y)
+z_z = circular_zero_conv(x, y)
+
+fig, ax = plt.subplots(3, 1, sharex=True, sharey=True)
+
+ax[0].grid()
+ax[1].grid()
+ax[2].grid()
+
+ax[0].stem(z_l.real)
+ax[1].stem(z_c.real)
+ax[2].stem(z_z.real)
+
+plt.savefig('./skotsugi/chapter03/q04.png')

skotsugi/chapter03/q04.py

@@ -0,0 +1,19 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+N = 10
+x = np.zeros(N)
+x[0] = 1
+
+i = np.arange(N)
+y = 0.2*x[i] + 0.2*x[i-1] + 0.2*x[i-2] + 0.2*x[i-3] + 0.2*x[i-4]
+
+fig, ax = plt.subplots(2, 1, sharex=True)
+
+ax[0].grid()
+ax[1].grid()
+
+ax[0].stem(x)
+ax[1].stem(y)
+
+plt.savefig('./skotsugi/chapter03/q05.png')

skotsugi/chapter03/q05.png

@@ -0,0 +1,24 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+N = 10
+x = np.zeros(N)
+x[0] = 1 # impulse
+
+def y(i: int):
+  if i < 0: return 0
+  return 0.3*y(i-1) + 0.4*x[i]
+
+i = np.arange(N)
+
+result = [y(i) for i in range(N)]
+
+fig, ax = plt.subplots(2, 1, sharex=True)
+
+ax[0].grid()
+ax[1].grid()
+
+ax[0].stem(x)
+ax[1].stem(result)
+
+plt.savefig('./skotsugi/chapter03/q06.png')

skotsugi/chapter03/q05.py

@@ -0,0 +1,45 @@
+import numpy as np
+
+def diff_equation_n(i, a, b, x):
+  N = len(a)
+  M = len(b)
+
+  if i < 0: 
+    return 0
+
+  y_sum = 0
+  x_sum = 0
+  for k in range(1, N):
+    j = i - k
+    y_sum += -a[k] * diff_equation_n(j, a, b, x)
+    
+  for k in range(M):
+    j = i - k
+    x_sum += (b[k] * x[j])
+  return y_sum + x_sum
+
+def diff_equation(a, b, x):
+  L = len(x)
+  return [diff_equation_n(i, a, b, x) for i in range(L)]
+
+
+if __name__ == "__main__":
+  import matplotlib.pyplot as plt
+
+  N = 10
+  x = np.zeros(N)
+  x[0] = 1 # impulse
+
+  a = [ 1, -0.3 ]
+  b = [ 0.4 ]
+  y = diff_equation(a, b, x)
+
+  fig, ax = plt.subplots(2, 1, sharex=True)
+
+  ax[0].grid()
+  ax[1].grid()
+
+  ax[0].stem(x)
+  ax[1].stem(y)
+
+  plt.savefig('./skotsugi/chapter03/q07.png')

skotsugi/chapter03/q06.png

@@ -0,0 +1,13 @@
+import numpy as np
+
+def freqz_w(w, a, b):
+  N = len(a)
+  M = len(b)
+
+  exp_a = np.array([ np.exp(-1j * w * i) for i in range(N)])
+  exp_b = np.array([ np.exp(-1j * w * i) for i in range(M)])
+
+  return np.sum(exp_b * b) / (np.sum(exp_a[1:N] * a[1:N]) + 1)
+
+def freqz(w, a, b):
+  return [freqz_w(i, a, b) for i in w]

skotsugi/chapter03/q06.py

@@ -0,0 +1,29 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q07 import diff_equation
+from q08 import freqz
+
+# impulse
+L = 10
+x = np.zeros(L)
+x[0] = 1
+a = [0]
+b = [0.2, 0.2, 0.2, 0.2, 0.2]
+
+y = diff_equation(a, b, x)
+
+fs = 16000
+N = 400
+f = np.arange(0, 1, 1 / N) * fs
+w = 2 * np.pi * f / fs
+
+result = freqz(w, a, b)
+
+A = np.abs(result)
+theta = np.angle(result)
+theta[A < 10**-5] = 0
+
+fig, ax = plt.subplots(1, 2)
+ax[0].stem(A)
+ax[1].stem(theta)
+plt.savefig('./skotsugi/chapter03/q09.png')

skotsugi/chapter03/q07.png

@@ -0,0 +1,28 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q07 import diff_equation
+from q08 import freqz
+
+N = 10
+x = np.zeros(N)
+x[0] = 1 # impulse
+
+a = [ 1, -0.3 ]
+b = [ 0.4 ]
+y = diff_equation(a, b, x)
+
+fs = 16000
+N = 400
+f = np.arange(0, 1, 1 / N) * fs
+w = 2 * np.pi * f / fs
+
+result = freqz(w, a, b)
+
+A = np.abs(result)
+theta = np.angle(result)
+theta[A < 10**-5] = 0
+
+fig, ax = plt.subplots(1, 2)
+ax[0].plot(A)
+ax[1].plot(theta)
+plt.savefig('./skotsugi/chapter03/q10.png')