Merge pull request #11 from woodnx/skotsugi/chapter02

Add chapter02

Author: taishi-n

skotsugi/chapter02/q01.py

@@ -0,0 +1,26 @@
+import numpy as np
+
+def dft(x):
+  N = len(x)
+  X = np.zeros(N, dtype = 'complex_')
+
+  for k in range(N):
+    for n in range(N):
+      X[k] += x[n] * np.exp(-1j * 2 * np.pi * k * n / N)
+  
+  return X
+
+def idft(X): 
+  N = len(X)
+  x = np.zeros(N, dtype = 'complex_')
+
+  for n in range(N):
+    for k in range(N):
+      x[k] += X[k] * np.exp(1j * 2 * np.pi * k * n / N) / N
+  
+  return x
+
+if __name__ == "__main__":
+  d = dft([1, 0, 0, 0])
+  D = idft(d)
+  print(d, D)

skotsugi/chapter02/q02.py

@@ -0,0 +1,5 @@
+from q01 import dft
+
+delta = [1, 0, 0, 0, 0, 0, 0, 0]
+
+print(dft(delta))

skotsugi/chapter02/q03.png

@@ -0,0 +1,22 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q01 import dft, idft
+
+delta = [1, 0, 0, 0, 0, 0, 0, 0]
+d = dft(delta)
+D = idft(d)
+
+print(D)
+
+fig, axs = plt.subplots(2)
+
+plt.ylim(0, 1)
+
+axs[0].stem(np.real(D))
+axs[0].set_title('Real Part')
+axs[1].stem(np.imag(D))
+axs[1].set_title('Imaginary Part')
+
+fig.tight_layout()
+
+plt.savefig('./skotsugi/chapter02/q03.png')

skotsugi/chapter02/q03.py

@@ -0,0 +1,26 @@
+import matplotlib.pyplot as plt
+
+def show_spectrum(filename, amp, phs): 
+  fig, axs = plt.subplots(1, 2)
+
+  axs[0].stem(amp)
+  axs[0].set_title('Amplitude')
+  axs[1].stem(phs)
+  axs[1].set_title('Phase')
+
+  fig.tight_layout()
+
+  plt.savefig(filename)
+
+if __name__ == "__main__":
+  import numpy as np
+  from q01 import dft
+
+  delta = [1, 0, 0, 0, 0, 0, 0, 0]
+  d = dft(delta)
+
+  amp = np.abs(d)
+  phs = np.angle(d)
+
+  show_spectrum('./skotsugi/chapter02/q04.png', amp, phs)
+  

skotsugi/chapter02/q04.png

@@ -0,0 +1,10 @@
+import numpy as np
+from q01 import dft
+
+delta = [1, 0, 0, 0, 0, 0, 0, 0]
+
+d = dft(delta)
+f = np.fft.fft(delta)
+
+diff = d - f
+print(diff)

skotsugi/chapter02/q04.py

@@ -0,0 +1,33 @@
+import numpy as np
+
+def sin(A: float, f: float, fs: float, sec: float):
+  n = np.arange(0, fs*sec) / fs
+  return A * np.sin(2 * np.pi * f * n)
+
+if __name__ == "__main__":
+  import matplotlib.pyplot as plt
+
+  fs = 16000
+  sec = 3.0
+  N = int(fs*sec)
+  
+  y = sin(1.0, 440, fs, sec)
+
+  Y = np.fft.fft(y)
+  freq = np.fft.fftfreq(N, d=1/fs)
+
+  amp = np.abs(Y)
+  phs = np.angle(Y)
+  phs[amp < 1e-5] = 0 # これをしないと，絶対値がすごく小さな値の偏角も計算されてしまい，でたらめな位相スペクトルが出てしまう
+  
+  fig, axs = plt.subplots(1, 2)
+
+  axs[0].stem(freq, amp)
+  axs[0].set_xlim([0, max(freq)])
+  axs[0].set_title('Amplitude')
+  axs[1].stem(freq, phs)
+  axs[1].set_xlim([0, max(freq)])
+  axs[1].set_title('Phase')
+
+  fig.tight_layout()
+  plt.savefig('./skotsugi/chapter02/q06.png')

skotsugi/chapter02/q05.py

@@ -0,0 +1,12 @@
+import numpy as np
+
+def hamming(N: int): 
+  return [ 0.54 - 0.46 * np.cos(2 * np.pi * n / (N - 1)) for n in range(int(N)) ]
+
+if __name__ == "__main__":
+  import matplotlib.pyplot as plt
+  
+  w = hamming(16000 * 3.0)
+  
+  plt.plot(w)
+  plt.savefig('./skotsugi/chapter02/q07.png')

skotsugi/chapter02/q06.png

@@ -0,0 +1,17 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q06 import sin
+from q07 import hamming
+
+# paramaters
+fs = 16000
+sec = 3.0
+f = 440
+N = fs * sec
+
+n = np.arange(0, N) / fs
+y = sin(1.0, f, fs, sec) * hamming(N)
+
+# plt.xlim(0, 0.05)
+plt.plot(n, y)
+plt.savefig('./skotsugi/chapter02/q08.png')

skotsugi/chapter02/q06.py

@@ -0,0 +1,28 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q07 import hamming
+
+fs = 16000
+sec = 3
+N = fs * sec
+wN = 2 ** 8
+
+w = hamming(N)
+W = np.fft.fft(w)
+
+amp = 20 * np.log10(np.abs(W))
+#amp = np.abs(W)
+phs = np.angle(W)
+phs[amp < 1e-5] = 0 
+
+freq = np.fft.fftfreq(int(N), d=1/fs)
+
+fig, axs = plt.subplots(1, 2)
+axs[0].stem(freq, amp)
+axs[0].set_title('Amplitude')
+axs[1].stem(freq, phs)
+axs[1].set_title('Phase')
+
+fig.tight_layout()
+
+plt.savefig('./skotsugi/chapter02/q09.png')

skotsugi/chapter02/q07.png

@@ -0,0 +1,37 @@
+import numpy as np
+
+def conv(X, Y):
+  N = len(X) 
+  Z = np.zeros(N, dtype = 'complex_')
+  k = np.arange(N)
+
+  for i in range(N):
+    Z[k] += X[i] * Y[k - i]
+
+  return Z
+
+
+if __name__ == "__main__":
+  import matplotlib.pyplot as plt
+  from q06 import sin
+  from q07 import hamming
+
+  fs = 16000
+  sec = 3
+  N = int(fs * sec)
+
+  x = sin(1, 440, fs, sec)
+  X = np.fft.fft(x)
+
+  y = hamming(N)
+  Y = np.fft.fft(y)
+
+  Z = conv(X, Y)
+  z = np.fft.ifft(Z) / N
+
+  print(z)
+
+  n = np.arange(0, N) / fs
+  #plt.xlim(0, 0.05)
+  plt.plot(n, z.real)
+  plt.savefig('./skotsugi/chapter02/q10.png')