Finished notch filter lab part 8

6th-Semester-Spring-2024/DSP/Labs/NotchFilter/dev_a.dat

@@ -7998,6 +7998,3 @@
   -6.3400000e+02
   -8.5800000e+02
   -1.0510000e+03
-
-
-

6th-Semester-Spring-2024/DSP/Labs/NotchFilter/notch_filter_8.py

@@ -0,0 +1,91 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.signal
+import sounddevice as sd
+
+
+f_s = 8000
+f_nyquist = f_s//2
+
+def normalize_signal(signal):
+    min_amp = np.min(signal)
+    normalized_signal = signal - min_amp
+    max_amp = np.max(normalized_signal)
+    normalized_signal *= 2/max_amp
+    normalized_signal -= 1
+    return normalized_signal
+
+def notch_filter(signal, f_0, f_s, r=0.99):
+    # Filter numerator coefficients
+    b_0 = 1
+    b_1 = -2*np.cos(2*np.pi*f_0/f_s)
+    b_2 = 1
+    b = [b_0, b_1, b_2]
+
+    # Filter denominator coefficients
+    a_0 = 1
+    a_1 = -2*r*np.cos(2*np.pi*f_0/f_s)
+    a_2 = r**2
+    a = [a_0, a_1, a_2]
+    
+    # Apply the filter
+    return scipy.signal.lfilter(b, a, signal)
+
+with open('dev_a.dat') as sound_data_file:
+    sound_data_strings = sound_data_file.readlines()
+    sound_data = []
+    for data_string in sound_data_strings:
+        sound_data.append(eval(data_string.strip()))
+    sound_data = np.array(sound_data)
+    
+    n = np.arange(len(sound_data))
+    omega = np.arange(f_nyquist)
+    
+
+    # Play original signal
+    #sd.play(normalize_signal(sound_data), samplerate=f_s, blocking=True)
+
+    # Plot signal frequency spectrum
+
+    freq_mag = 20*np.log10(np.abs(np.fft.fft(sound_data)))
+    plt.figure(figsize=(16,9))
+    plt.scatter(omega,freq_mag[:f_nyquist])
+    plt.xlabel("Frequency [Hz]")
+    plt.ylabel("Magnitude [dB]")
+    plt.title("Frequency Spectrum without Noise")
+    plt.savefig("dev_a_freq_no_noise.png")
+
+    # Calculate noise
+    noise_amp = np.max(np.abs(sound_data))/2
+    f_noise = 1000
+    noise = noise_amp * np.cos(2*np.pi*f_noise*n/f_s)
+
+    # Add noise
+    noisy_sound_data = noise+sound_data
+
+    # Play signal+noise
+    #sd.play(normalize_signal(noisy_sound_data), samplerate=f_s, blocking=True)
+
+    # Plot signal+noise frequency spectrum
+    freq_mag = 20*np.log10(np.abs(np.fft.fft(noisy_sound_data)))
+    plt.figure(figsize=(16,9))
+    plt.scatter(omega,freq_mag[:f_nyquist])
+    plt.xlabel("Frequency [Hz]")
+    plt.ylabel("Magnitude [dB]")
+    plt.title("Frequency Spectrum with Noise")
+    plt.savefig("dev_a_freq_noise.png")
+
+    # Apply filter
+    filtered_sound_data = notch_filter(noisy_sound_data, f_0=f_noise, f_s=f_s)
+
+    # Play filtered signal+noise
+    #sd.play(normalize_signal(filtered_sound_data), samplerate=f_s, blocking=True)
+
+    # Plot filtered signal+noise
+    freq_mag = 20*np.log10(np.abs(np.fft.fft(filtered_sound_data)))
+    plt.figure(figsize=(16,9))
+    plt.scatter(omega,freq_mag[:f_nyquist])
+    plt.xlabel("Frequency [Hz]")
+    plt.ylabel("Magnitude [dB]")
+    plt.title("Frequency Spectrum with Noise Filtered")
+    plt.savefig("dev_a_freq_filtered.png")