２章を追加 q06,q10についてはうまくいっていないため修正予定

Author: kuu8902

kkoiso/chapter02/q01.py

@@ -0,0 +1,24 @@
+import numpy as np
+import matplotlib.pyplot as plt
+def DFT(x):
+    N = len(x)
+    X = [0] * N
+    for k in range(N):
+        s=0
+        for n in range(N): 
+            a = ((-2j) * np.pi * k * n) / N
+            s = s + x[n] * np.exp(a)
+           
+        X[k] = s
+    return X
+
+def IDFT(X):
+    N = len(X)
+    x = [0] * N
+    for n in range(N):
+        s=0
+        for k in range(N):
+            angle = (2j * np.pi * k * n) / N
+            s = s + X[k] * np.exp(angle)
+        x[n] =s/N
+    return x

kkoiso/chapter02/q02.py

@@ -0,0 +1,10 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q01 import DFT,IDFT
+
+x = np.array([1, 0, 0, 0, 0, 0, 0, 0])
+X = DFT(x)
+
+
+print(X)
+

kkoiso/chapter02/q03.py

@@ -0,0 +1,28 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q01 import DFT,IDFT
+
+x = np.array([1, 0, 0, 0, 0, 0, 0, 0])
+X = DFT(x)
+x_idft = IDFT(X)
+print(x_idft)
+idft_real = [val.real for val in x_idft]#実部を取りだす
+idft_imag = [val.imag for val in x_idft]#虚部を取り出す
+
+#実部をプロット
+plt.figure(figsize=(12, 6))
+plt.subplot(2, 1, 1)
+plt.stem(idft_real)
+plt.title('IDFT(Real Part)')
+plt.xlabel('Sample')
+plt.ylabel('Amplitude')
+
+
+plt.subplot(2, 1, 2)
+plt.stem(idft_imag)
+plt.title('IDFT Result (Imaginary Part)')
+plt.xlabel('Sample')
+plt.ylabel('Amplitude')
+
+plt.tight_layout()
+plt.show()

kkoiso/chapter02/q04.py

@@ -0,0 +1,25 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q01 import DFT,IDFT
+
+x = np.array([1, 0, 0, 0, 0, 0, 0, 0])
+X = DFT(x)
+title = "DFT impulse"
+
+plt.figure(figsize=(10, 4))
+plt.subplot(2, 1, 1)
+plt.plot(np.abs(X))
+plt.title(title + ' Amplitude Spectrum')
+plt.xlabel('Frequency')
+plt.ylabel('Amplitude')
+plt.grid(True)
+    
+plt.subplot(2, 1, 2)
+plt.plot(np.angle(X))
+plt.title(title + ' Phase Spectrum')
+plt.xlabel('Frequency')
+plt.ylabel('Phase')
+plt.grid(True)
+    
+plt.tight_layout()
+plt.show()

kkoiso/chapter02/q05.py

@@ -0,0 +1,11 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q01 import DFT,IDFT
+
+x = np.array([1, 0, 0, 0, 0, 0, 0, 0])
+X = DFT(x)
+
+
+np_X = np.fft.fft(x)
+print("Custom DFT:", X)
+print("Np FFT:", np_X)

kkoiso/chapter02/q06.py

@@ -0,0 +1,45 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+a = 1.0
+freq= 440 
+samp_rate = 16000 
+d = 3  
+
+t = np.linspace(0, d, int(samp_rate * d), endpoint=False)
+
+
+sin_wave = a * np.sin(2 * np.pi * freq * t)
+
+X = np.fft.fft(sin_wave)
+
+#振幅をデシベル表記
+amplitude_spectrum = 20 * np.log10(np.abs(X))
+
+# 位相スペクトルの計算
+phase_spectrum = np.angle(X)
+
+# 対応する周波数を計算
+frequencies = np.fft.fftfreq(len(X), 1/samp_rate)
+
+
+plt.figure(figsize=(12, 6))
+# 振幅スペクトルのプロット
+plt.subplot(2, 1, 1)
+plt.plot(frequencies,amplitude_spectrum)
+plt.title('Amplitude Spectrum')
+plt.xlabel('Frequency (Hz)')
+plt.ylabel('Amplitude (dB)')
+plt.grid(True)
+
+# 位相スペクトルのプロット
+plt.subplot(2, 1, 2)
+plt.plot(frequencies, phase_spectrum)
+plt.title('Phase Spectrum')
+plt.xlabel('Frequency (Hz)')
+plt.ylabel('Phase (radians)')
+plt.grid(True)
+
+plt.tight_layout()
+plt.show()

kkoiso/chapter02/q07.py

@@ -0,0 +1,6 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def Hamming_Window(N):
+    n = np.arange(N)
+    return 0.54 - 0.46 * np.cos((2 * np.pi * n )/ (N - 1))

kkoiso/chapter02/q08.py

@@ -0,0 +1,25 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q07 import Hamming_Window
+#正弦波を生成
+a = 1.0
+freq= 440 
+samp_rate = 16000 
+d = 3  
+
+t = np.linspace(0, d, int(samp_rate * d), endpoint=False)
+
+
+sin_wave = a * np.sin(2 * np.pi * freq * t)
+#ハミング窓を生成し、正弦波にかける
+signal_length = len(sin_wave)
+hamming_window_signal = sin_wave * Hamming_Window(signal_length)
+
+#プロット
+plt.figure(figsize=(10, 4))
+plt.plot(hamming_window_signal)
+plt.title('Signal with Hamming Window')
+plt.xlabel('Sample')
+plt.ylabel('Amplitude')
+plt.grid(True)
+plt.show()

kkoiso/chapter02/q09.py

@@ -0,0 +1,19 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q07 import Hamming_Window
+#正弦波を生成
+a = 1.0
+freq= 440 
+samp_rate = 16000 
+d = 3  
+
+t = np.linspace(0, d, int(samp_rate * d), endpoint=False)
+
+
+sin_wave = a * np.sin(2 * np.pi * freq * t)
+
+signal_length = len(sin_wave)
+#正弦波と同じサイズのハミング窓を生成し、DFT
+hamming_window_spectrum = np.fft.fft(Hamming_Window(signal_length))
+
+print(hamming_window_spectrum)

kkoiso/chapter02/q10.py

@@ -0,0 +1,41 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from q07 import Hamming_Window
+from q01 import DFT,IDFT
+#正弦波を生成
+a = 1.0
+freq= 440 
+samp_rate = 16000 
+d = 3  
+
+t = np.linspace(0, d, int(samp_rate * d), endpoint=False)
+
+
+sin_wave = a * np.sin(2 * np.pi * freq * t)
+#正弦波をDFT
+sin_wave_DFT = np.fft.fft(sin_wave)
+
+
+#正弦波と同じサイズのハミング窓を生成し、DFT
+signal_length = len(sin_wave)
+hamming_window_spectrum = np.fft.fft(Hamming_Window(signal_length))
+
+#巡回畳み込み
+Z=[0] * signal_length
+for k in range(signal_length):
+    s = 0
+    for n in range(signal_length):
+        s =s +sin_wave_DFT[n]*hamming_window_spectrum[k-n]
+    Z[k] = s
+z_idft = np.fft.ifft(Z)
+print(z_idft)
+
+
+#プロット
+plt.figure(figsize=(10, 4))
+plt.plot(z_idft)
+plt.title('Signal with Hamming Window')
+plt.xlabel('Sample')
+plt.ylabel('Amplitude')
+plt.grid(True)
+plt.show()