Merge pull request #12 from 153hashimoto/khashimoto/chapter02

Add chapter02

khashimoto/chapter02/q01.py

@@ -0,0 +1,28 @@
+# DFT/IDFT の実装: N 点の信号の離散フーリエ変換 (DFT)とその逆変換(IDFT)を計算する関数を実装せよ．
+
+import numpy as np
+
+def dft(x):
+    N = len(x)
+    dft_X = np.zeros(N, dtype=complex)
+    for k in range(N):
+        for n in range(N):
+            dft_X[k] += x[n] * np.exp(-1j * 2 * np.pi * k * n / N)
+    return dft_X
+
+def idft(X):
+    N = len(X)
+    idft_x = np.zeros(N, dtype=complex)
+    for n in range(N):
+        for k in range(N):
+            idft_x[n] += X[k] * np.exp(1j * 2 * np.pi * k * n / N)
+        idft_x[n] /= N
+    return idft_x
+
+
+# --------------- 確認----------------------
+x_1 = np.array([1, 0, 0, 0, 0])
+X_1 = dft(x_1)
+print(X_1)
+x_2 = idft(x_1)
+print(x_2)

khashimoto/chapter02/q02.py

@@ -0,0 +1,15 @@
+# DFT の確認: 1.で実装した関数を用いて 8 点の単位インパルス信号 の DFT を計算せよ.
+
+import numpy as np
+
+def dft(x):
+    N = len(x)
+    dft_X = np.zeros(N, dtype=complex)
+    for k in range(N):
+        for n in range(N):
+            dft_X[k] += x[n] * np.exp(-1j * 2 * np.pi * k * n / N)
+    return dft_X
+
+delta = np.array([1, 0, 0, 0, 0, 0, 0, 0])
+X = dft(delta)
+print(X)

khashimoto/chapter02/q03.py

@@ -0,0 +1,56 @@
+# IDFT の確認: 2.の結果の IDFT を計算しプロットせよ.
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+def dft(x):
+    N = len(x)
+    dft_X = np.zeros(N, dtype=complex)
+    for k in range(N):
+        for n in range(N):
+            dft_X[k] += x[n] * np.exp(-1j * 2 * np.pi * k * n / N)
+    return dft_X
+
+def idft(X):
+    N = len(X)
+    idft_x = np.zeros(N, dtype=complex)
+    for n in range(N):
+        for k in range(N):
+            idft_x[n] += X[k] * np.exp(1j * 2 * np.pi * k * n / N)
+        idft_x[n] /= N
+    return idft_x
+
+
+
+delta = np.array([1, 0, 0, 0, 0, 0, 0, 0])
+X = dft(delta)    # 2. の結果
+idft_x = idft(X)
+plt.stem(idft_x.real)
+plt.show()
+
+
+# ---------------- 確認 ----------------------
+'''
+x = np.arange(0, 2 * np.pi, 0.1)
+
+fig = plt.figure()
+plt.subplot(5,1,1)
+y1 = np.sin(x)
+plt.plot(y1)
+
+plt.subplot(5,1,2)
+Y = dft(y1)
+plt.plot(Y.real)
+
+plt.subplot(5,1,3)
+plt.plot(Y.imag)
+
+y2 = idft(Y)
+plt.subplot(5,1,4)
+plt.plot(y2.real)
+
+plt.subplot(5,1,5)
+plt.plot(y2.imag)
+
+plt.show()
+'''

khashimoto/chapter02/q04.py

@@ -0,0 +1,22 @@
+# スペクトルの確認: 2.の結果の振幅スペクトルおよび位相スペクトル をプロットせよ．
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+def dft(x):
+    N = len(x)
+    dft_X = np.zeros(N, dtype=complex)
+    for k in range(N):
+        for n in range(N):
+            dft_X[k] += x[n] * np.exp(-1j * 2 * np.pi * k * n / N)
+    return dft_X
+
+delta = np.array([1, 0, 0, 0, 0, 0, 0, 0])
+X = dft(delta)    # 2. の結果
+
+plt.subplot(2,1,1)
+plt.stem(abs(X))
+
+plt.subplot(2,1,2)
+plt.stem(np.angle(X))
+plt.show()

khashimoto/chapter02/q05.py

@@ -0,0 +1,8 @@
+# 既存の実装との比較: 8 点の単位インパルス信号の DFT を numpy.fft.fft 関数を用いて計算し，2.の結果と比較せよ．なお，これ以降 DFT を計算する場合は numpy.fft.fft 関数を使用してよい.
+
+import numpy as np
+
+delta = np.array([1, 0, 0, 0, 0, 0, 0, 0])
+X = np.fft.fft(delta)
+print(X)
+# 2. と同じ結果となった

khashimoto/chapter02/q06.py

@@ -0,0 +1,23 @@
+#  正弦波のスペクトル: 第 1 章 1.で作成した信号の DFT を計算し，振幅スペクトルと位相スペクトルをプロットせよ．（プロットする際はデシベル表記にするために 20 * log10(np.abs(X)) などのようにしてデシベル表記になおしてから計算すると見やすいです．）
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+A = 1    # 振幅
+f = 440    # 周波数
+fs = 16000    # サンプリング周波数
+T = 3
+t = np.arange(fs * T) / fs
+x = A * np.sin(2 * np.pi * f * t)    # 第１章１. で作成した信号
+
+X = np.fft.fft(x)
+
+
+# 振幅スペクトルと位相スペクトルをプロット
+plt.subplot(2,1,1)
+plt.plot(20 * np.log10(abs(X)))
+#plt.stem(abs(X))
+
+plt.subplot(2,1,2)
+plt.plot(np.angle(X))
+plt.show()

khashimoto/chapter02/q07.py

@@ -0,0 +1,14 @@
+# 窓関数: 次式で定義される 点の Hamming 窓を作成する関数を実装せよ.
+
+import numpy as np
+
+def hamming(N):
+    w = np.zeros(N)
+    for n in range(N):
+        w[n] = 0.54-0.46 * np.cos(2 * np.pi * n / (N-1))
+    return w
+
+
+# ------------- 確認 ----------------
+w = hamming(5)
+print(w)

khashimoto/chapter02/q08.py

@@ -0,0 +1,24 @@
+# 窓関数の確認: 第 1 章 1.で作成した信号に対して 6.の窓関数を適用した信号をプロットせよ．
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+def hamming(N):
+    w = np.zeros(N)
+    for n in range(N):
+        w[n] = 0.54-0.46 * np.cos(2 * np.pi * n / (N-1))
+    return w
+
+
+A = 1    # 振幅
+f = 440    # 周波数
+fs = 16000    # サンプリング周波数
+T = 3
+t = np.arange(fs * T) / fs
+x = A * np.sin(2 * np.pi * f * t)    # 第１章１. で作成した信号
+
+w = hamming(fs*T)
+x2 = x * w
+
+plt.plot(t, x2)
+plt.show()

khashimoto/chapter02/q09.py

@@ -0,0 +1,16 @@
+# 窓関数の特性: 第 1 章 1.で作成した信号と同じ長さの Hamming 窓を作成し，DFT を計算せよ．
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+def hamming(N):
+    w = np.zeros(N)
+    for n in range(N):
+        w[n] = 0.54-0.46 * np.cos(2 * np.pi * n / (N-1))
+    return w
+
+fs = 16000    # サンプリング周波数
+T = 3    # 時間
+w = hamming(fs*T)
+w = np.fft.fft(w)
+print(w)

khashimoto/chapter02/q10.py

@@ -0,0 +1,33 @@
+# 窓関数とスペクトルの関係: 
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+def hamming(N):
+    w = np.zeros(N)
+    for n in range(N):
+        w[n] = 0.54-0.46 * np.cos(2 * np.pi * n / (N-1))
+    return w
+
+
+A = 1    # 振幅
+f = 440    # 周波数
+fs = 16000    # サンプリング周波数
+T = 3
+t = np.arange(fs * T) / fs
+x = A * np.sin(2 * np.pi * f * t)    # 第１章１. で作成した信号
+X = np.fft.fft(x)
+
+w = hamming(fs*T)
+Y = np.fft.fft(w)    # ９. の結果
+
+# X[k]とY[k]の巡回畳み込み
+Z = np.zeros(fs*T, dtype=complex)
+for k in range(fs*T):
+    for n in range(fs*T):
+        Z[k] += X[n] * Y[k-n]
+
+idft_z = np.fft.ifft(Z) / (fs * T)
+
+plt.plot(t, idft_z.real)
+plt.show()