Merge https://github.com/Yuto-Yamamoto01/asp-tutorial-2024 into yyamamoto/chapter02

Author: Yuto-Yamamoto01

khashimoto/chapter03/q01.py

@@ -0,0 +1,15 @@
+# 線形畳み込み
+
+import numpy as np
+
+
+def LinearConv(x, h):
+    N = x.size
+    z = np.zeros(2 * N - 1)
+    for n in range(2 * N - 1):
+        for k in range(N):
+            if (n - k) < 0 or (n - k) > (N - 1):
+                z[n] += 0
+            else:
+                z[n] += x[k] * h[n - k]
+    return z

khashimoto/chapter03/q02.py

@@ -0,0 +1,11 @@
+# 巡回畳み込み
+
+import numpy as np
+
+
+def CircularConv(x, h):
+    N = x.size
+    z = np.zeros(N)
+    for k in range(N):
+        z += x[k] * h[(np.arange(N) - k) % N]
+    return z

khashimoto/chapter03/q03.py

@@ -0,0 +1,14 @@
+# 零詰め+巡回畳み込み
+
+import numpy as np
+
+
+def ZeroPaddingCircularConv(x, h):
+    N1 = x.size
+    x = np.hstack((x, np.zeros(N1 - 1)))
+    h = np.hstack((h, np.zeros(N1 - 1)))
+    N = 2 * N1 - 1
+    z = np.zeros(N)
+    for k in range(N):
+        z += x[k] * h[(np.arange(N) - k) % N]
+    return z

khashimoto/chapter03/q04.py

@@ -0,0 +1,15 @@
+# 各種畳み込みの関係
+
+import numpy as np
+from q01 import LinearConv
+from q02 import CircularConv
+from q03 import ZeroPaddingCircularConv
+
+x = np.array([4, 3, 2, 1])
+h = np.array([1, 0, -1, 0])
+print("線形畳み込み:")
+print(LinearConv(x, h))
+print("巡回畳み込み:")
+print(CircularConv(x, h))
+print("零詰めを行った巡回畳み込み:")
+print(ZeroPaddingCircularConv(x, h))

khashimoto/chapter03/q05.py

@@ -0,0 +1,12 @@
+# 差分方程式（再帰なし）
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+x = np.array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0])
+y = np.zeros(10)
+n = np.arange(10)
+y[n] = 0.2 * x[n] + 0.2 * x[n - 1] + 0.2 * x[n - 2] + 0.2 * x[n - 3] + 0.2 * x[n - 4]
+
+plt.stem(y)
+plt.show()

khashimoto/chapter03/q06.py

@@ -0,0 +1,12 @@
+# 差分方程式（再帰あり）
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+x = np.array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0])
+y = np.zeros(10)
+for n in range(10):
+    y[n] = 0.3 * y[n - 1] + 0.4 * x[n]
+
+plt.stem(y)
+plt.show()

khashimoto/chapter03/q07.py

@@ -0,0 +1,17 @@
+# 差分方程式（一般系）
+
+import numpy as np
+
+
+def DifferenceEquation(a, b, x):
+    N = a.size - 1
+    M = b.size - 1
+    L = x.size
+
+    y = np.zeros(L)
+    ka = np.arange(1, N + 1)
+    kb = np.arange(M + 1)
+    for n in range(L):
+        y[n] = (-sum(a[ka] * y[n - ka]) + sum(b[kb] * x[n - kb])) / a[0]
+
+    return y

khashimoto/chapter03/q08.py

@@ -0,0 +1,15 @@
+# 周波数応答
+
+import numpy as np
+
+
+def FreqRes(a, b, fs, f):
+    N = a.size - 1
+    M = b.size - 1
+
+    ka = np.arange(1, N + 1)
+    kb = np.arange(M + 1)
+    H = sum(b[kb] * np.exp(-1j * 2 * np.pi * f / fs * kb)) / (
+        1 + sum(a[ka] * np.exp(-1j * 2 * np.pi * f / fs * ka))
+    )
+    return H

khashimoto/chapter03/q09.py

@@ -0,0 +1,22 @@
+# 周波数応答の確認（再帰なし）
+
+import numpy as np
+import matplotlib.pyplot as plt
+from q08 import FreqRes
+
+fs = 16000
+N = 10
+f = np.arange(0, fs, 1 / N)
+
+a = np.array([0])
+b = np.array([0.2, 0.2, 0.2, 0.2, 0.2])
+H = np.zeros(f.size, dtype=complex)
+for i in range(f.size):
+    H[i] = FreqRes(a, b, fs, f[i])
+
+plt.subplot(2, 1, 1)
+plt.plot(abs(H))
+
+plt.subplot(2, 1, 2)
+plt.plot(np.angle(H))
+plt.show()

khashimoto/chapter03/q10.py

@@ -0,0 +1,22 @@
+# 周波数応答の確認（再帰あり）
+
+import numpy as np
+import matplotlib.pyplot as plt
+from q08 import FreqRes
+
+fs = 16000
+N = 10
+f = np.arange(0, fs, 1 / N)
+
+a = np.array([0.3])
+b = np.array([0.4])
+H = np.zeros(f.size, dtype=complex)
+for i in range(f.size):
+    H[i] = FreqRes(a, b, fs, f[i])
+
+plt.subplot(2, 1, 1)
+plt.plot(abs(H))
+
+plt.subplot(2, 1, 2)
+plt.plot(np.angle(H))
+plt.show()

kkoiso/chapter03/q01.py

@@ -0,0 +1,12 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def linear_conv(x, h):
+    N = len(x)
+    z = np.zeros(2 * N - 1)
+    for n in range(2 * N - 1):
+        for k in range(N):
+            if 0 <= n - k <= N - 1:
+                z[n] = z[n] + x[k] * h[n - k]
+    return z

kkoiso/chapter03/q02.py

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

kkoiso/chapter03/q03.py

@@ -0,0 +1,11 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q02 import circular_conv
+
+
+def zero_padded_circular_conv(x, h):
+    N = len(x)
+    x_padded = np.concatenate([x, np.zeros(N)])
+    h_padded = np.concatenate([h, np.zeros(N)])
+    z = circular_conv(x_padded, h_padded)
+    return z

kkoiso/chapter03/q04.py

@@ -0,0 +1,33 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q01 import linear_conv
+from q02 import circular_conv
+from q03 import zero_padded_circular_conv
+
+x = [4, 3, 2, 1]
+h = [1, 0, -1, 0]
+x_padded = np.concatenate([x, np.zeros(4)])
+
+z_linear = linear_conv(x, h)
+z_circular = circular_conv(x, h)
+z_zero_padded = zero_padded_circular_conv(x, h)
+
+plt.figure(figsize=(12, 8))
+
+plt.subplot(3, 1, 1)
+plt.stem(z_linear)
+plt.title("Linear Convolution")
+plt.grid(True)
+
+plt.subplot(3, 1, 2)
+plt.stem(z_circular)
+plt.title("Circular Convolution")
+plt.grid(True)
+
+plt.subplot(3, 1, 3)
+plt.stem(z_zero_padded)
+plt.title("Zero-padded Circular Convolution")
+plt.grid(True)
+
+plt.tight_layout()
+plt.show()

kkoiso/chapter03/q05.py

@@ -0,0 +1,27 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def difference_equations_no_recursion(x):
+    y = np.zeros(len(x))
+    for n in range(len(x)):
+        y[n] = 0.2 * x[n]
+        if n >= 1:
+            y[n] = y[n] + 0.2 * x[n - 1]
+        if n >= 2:
+            y[n] = y[n] + 0.2 * x[n - 2]
+        if n >= 3:
+            y[n] = y[n] + 0.2 * x[n - 3]
+        if n >= 4:
+            y[n] = y[n] + 0.2 * x[n - 4]
+    return y
+
+
+x = np.zeros(10)
+x[0] = 1
+y = difference_equations_no_recursion(x)
+
+plt.stem(y)
+plt.title("Difference Equation (Non-recursive)")
+plt.grid(True)
+plt.show()

kkoiso/chapter03/q06.py

@@ -0,0 +1,22 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def difference_equations_recursion(x):
+    y = np.zeros(len(x))
+    for n in range(len(x)):
+        if n == 0:
+            y[n] = 0.4 * x[n]
+        else:
+            y[n] = 0.3 * y[n - 1] + 0.4 * x[n]
+    return y
+
+
+x = np.zeros(10)
+x[0] = 1
+y = difference_equations_recursion(x)
+
+plt.stem(y)
+plt.title("Difference Equation (Recursive)")
+plt.grid(True)
+plt.show()

kkoiso/chapter03/q07.py

@@ -0,0 +1,18 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def general_difference_equations(a, b, x):
+    N = len(a)
+    M = len(b)
+    L = len(x)
+    y = np.zeros(L)
+    for n in range(1, L + 1):
+        for k in range(1, N + 1):
+            if n >= k:
+                y[n] = y[n] - a[k] * y[n - k]
+        for k in range(M + 1):
+            if n >= k:
+                y[n] = y[n] + b[k] * x[n - k]
+        y[n] = y[n] / a[0]
+    return y

kkoiso/chapter03/q08.py

@@ -0,0 +1,20 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def frequency_response(a, b, fs):
+    omega = np.linspace(0, 2 * np.pi, fs, endpoint=False)
+
+    num = np.zeros(len(omega), dtype=np.complex128)
+    den = np.zeros(len(omega), dtype=np.complex128)
+
+    for k in range(len(b)):
+        num = num + b[k] * np.exp(-1j * (omega) * k)
+
+    den = den + 1
+    for k in range(1, len(a)):
+        den = den + a[k] * np.exp(-1j * (omega) * k)
+
+    H = num / den
+
+    return omega, H

kkoiso/chapter03/q09.py

@@ -0,0 +1,32 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q08 import frequency_response
+
+a = [1]
+b = [0.2, 0.2, 0.2, 0.2, 0.2]
+# b = [0.33, 0.33, 0.33]
+fs = 16000
+N = 16000
+w, H = frequency_response(a, b, fs)
+w = w[:N]
+H = H[:N]
+
+
+plt.figure(figsize=(12, 6))
+
+plt.subplot(2, 1, 1)
+plt.plot(w / (2 * np.pi) * fs, np.abs(H))
+plt.title("Frequency Response - No-recursion (Amplitude)")
+plt.xlabel("Frequency (Hz)")
+plt.ylabel("Amplitude (dB)")
+plt.grid(True)
+
+plt.subplot(2, 1, 2)
+plt.plot(w / (2 * np.pi) * fs, np.angle(H))
+plt.title("Frequency Response - No-recursion (Phase)")
+plt.xlabel("Frequency (Hz)")
+plt.ylabel("Phase (radians)")
+plt.grid(True)
+
+plt.tight_layout()
+plt.show()

kkoiso/chapter03/q10.py

@@ -0,0 +1,29 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q08 import frequency_response
+
+a = [1, -0.3]
+b = [0.4]
+fs = 16000
+N = 16000
+w, H = frequency_response(a, b, fs)
+w = w[:N]
+H = H[:N]
+plt.figure(figsize=(12, 6))
+
+plt.subplot(2, 1, 1)
+plt.plot(w / (2 * np.pi) * fs, np.abs(H))
+plt.title("Frequency Response - Recursive (Ampliitude)")
+plt.xlabel("Frequency (Hz)")
+plt.ylabel("Amplitude (dB)")
+plt.grid(True)
+
+plt.subplot(2, 1, 2)
+plt.plot(w / (2 * np.pi) * fs, np.angle(H))
+plt.title("Frequency Response - Recursive (Phase)")
+plt.xlabel("Frequency (Hz)")
+plt.ylabel("Phase (radians)")
+plt.grid(True)
+
+plt.tight_layout()
+plt.show()