|
|
@@ -0,0 +1,109 @@
|
|
|
+import numpy as np
|
|
|
+import matplotlib.pyplot as plt
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+class CICFilter:
|
|
|
+ def __init__(self, decimation_factor, delay, order):
|
|
|
+ self.decimation_factor = decimation_factor
|
|
|
+ self.delay = delay
|
|
|
+ self.order = order
|
|
|
+
|
|
|
+ def integrate(self, input_signal):
|
|
|
+ output = np.copy(input_signal)
|
|
|
+ for _ in range(self.order):
|
|
|
+ output = np.cumsum(output)
|
|
|
+ return output
|
|
|
+
|
|
|
+ def comb(self, decimated_signal):
|
|
|
+ output = np.copy(decimated_signal)
|
|
|
+ for _ in range(self.order):
|
|
|
+ delayed = np.zeros_like(output)
|
|
|
+ delayed[self.delay:] = output[:-self.delay]
|
|
|
+ output = output - delayed
|
|
|
+ return output
|
|
|
+
|
|
|
+ def decimate(self, integrated_signal):
|
|
|
+ # Децимация: выбираем каждый N-й элемент
|
|
|
+ return integrated_signal[::self.decimation_factor]
|
|
|
+
|
|
|
+
|
|
|
+ def filter(self, input_signal):
|
|
|
+ # 1. Каскад интеграторов
|
|
|
+ integrated_signal = self.integrate(input_signal)
|
|
|
+
|
|
|
+ # 2. Децимация
|
|
|
+ decimated_signal = self.decimate(integrated_signal)
|
|
|
+
|
|
|
+ # 3. Каскад comb-фильтров
|
|
|
+ output_signal = self.comb(decimated_signal)
|
|
|
+
|
|
|
+ return output_signal
|
|
|
+
|
|
|
+ def get_filter_freq_resp(self, num_points=1024):
|
|
|
+ # Подаём дельта-функцию на вход фильтру
|
|
|
+ impulse = np.zeros(num_points)
|
|
|
+ impulse[0] = 1.0 # Дельта-функция
|
|
|
+ output_signal = self.filter(impulse)
|
|
|
+ # Вычисляем частотную характеристику с помощью БПФ
|
|
|
+ freq_response = np.fft.fft(output_signal, num_points)
|
|
|
+ # Нормализуем частотную характеристику
|
|
|
+ freq_response /= np.max(np.abs(freq_response))
|
|
|
+ # Создаем массив частот
|
|
|
+ freq = np.fft.fftfreq(num_points, d=1/(50e6/self.decimation_factor))
|
|
|
+ # Возвращаем частоты и амплитуду частотной характеристики
|
|
|
+ return freq, np.abs(freq_response)
|
|
|
+
|
|
|
+ # # Calculate coefficients for the FIR filter to compensate the CIC frequency response
|
|
|
+ # def get_compensation_coefficients(self, cic_freq_response, target_freq_response):
|
|
|
+
|
|
|
+
|
|
|
+def main():
|
|
|
+ decimation_factor = 4
|
|
|
+ delay = 8
|
|
|
+ order = 6
|
|
|
+ num_samples = 100000
|
|
|
+ # Input signal - multiplication of a 3 MHz sine and 6 MHz sine
|
|
|
+ fs = 50e6 # Sampling frequency
|
|
|
+ f1 = 7e6 # Frequency of first sine wave
|
|
|
+ f2 = 3e6 # Frequency of second sine wave
|
|
|
+ t = np.arange(num_samples) / fs
|
|
|
+ input_signal = np.sin(2 * np.pi * f1 * t) * np.sin(2 * np.pi * f2 * t)
|
|
|
+ # Plot the input signal in the frequency domain
|
|
|
+ plt.figure(figsize=(10, 6))
|
|
|
+ plt.magnitude_spectrum(input_signal, Fs=fs, scale='dB', color='C0')
|
|
|
+ plt.title('Input Signal in Frequency Domain')
|
|
|
+ plt.xlabel('Frequency [Hz]')
|
|
|
+ plt.ylabel('Magnitude [dB]')
|
|
|
+ plt.grid()
|
|
|
+ plt.tight_layout()
|
|
|
+ plt.show()
|
|
|
+ # Create CIC filter instance
|
|
|
+ cic_filter = CICFilter(decimation_factor, delay, order)
|
|
|
+ # Filter the input signal
|
|
|
+ output_signal = cic_filter.filter(input_signal)
|
|
|
+ # Plot the output signal in the frequency domain
|
|
|
+ plt.magnitude_spectrum(output_signal, Fs=fs/decimation_factor, scale='dB', color='C1')
|
|
|
+ plt.title('Output Signal in Frequency Domain')
|
|
|
+ plt.xlabel('Frequency [Hz]')
|
|
|
+ plt.ylabel('Magnitude [dB]')
|
|
|
+ plt.grid()
|
|
|
+ plt.tight_layout()
|
|
|
+ plt.show()
|
|
|
+ # Get frequency response of the CIC filter
|
|
|
+ freq, freq_response = cic_filter.get_filter_freq_resp()
|
|
|
+ # Plot the frequency response
|
|
|
+ plt.figure(figsize=(10, 6))
|
|
|
+ # Normalize frequency axis to 1
|
|
|
+ freq_normalized = freq[:len(freq)//2]
|
|
|
+ plt.plot(freq_normalized, 20 * np.log10(freq_response[:len(freq)//2]), label='CIC Filter Frequency Response (Normalized)')
|
|
|
+ plt.title('CIC Filter Frequency Response')
|
|
|
+ plt.xlabel('Frequency [Hz]')
|
|
|
+ plt.ylabel('Magnitude [dB]')
|
|
|
+ plt.grid()
|
|
|
+ plt.legend()
|
|
|
+ plt.tight_layout()
|
|
|
+ plt.show()
|
|
|
+
|
|
|
+if __name__ == "__main__":
|
|
|
+ main()
|