Created
May 23, 2021 15:30
-
-
Save z0gSh1u/4ee358abed59ac31461d5779ffad7882 to your computer and use it in GitHub Desktop.
FIR Filter
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # coding: utf-8 | |
| # 《数字信号处理》课程实验 | |
| # FIR数字滤波器设计 | |
| # 09017227 卓旭 | |
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| plt.rcParams['font.sans-serif'] = ['KaiTi'] # 指定默认字体 | |
| plt.rcParams['axes.unicode_minus'] = False | |
| import math | |
| import random | |
| from fft import * | |
| W_C = 0.2 * math.pi | |
| N = 32 | |
| ''' | |
| 未加窗的无限长h_d(n) | |
| ''' | |
| def h_d(w_c): | |
| alpha = (N - 1) / 2. # 系统群时延 | |
| vals = [] | |
| for i in range(N * 3): | |
| vals.append(math.sin(w_c * (i - alpha)) / float(math.pi * (i - alpha))) | |
| return vals | |
| ''' | |
| 矩形窗 | |
| ''' | |
| def rect_window(one_len): | |
| res = [] | |
| for i in range(one_len): | |
| res.append(1.) | |
| return res | |
| ''' | |
| Hamming窗 | |
| ''' | |
| def Hamming_window(length): | |
| res = [] | |
| for i in range(length): | |
| res.append(0.54 - 0.46 * math.cos(2 * math.pi * i / (length - 1))) | |
| return res | |
| ''' | |
| Barlett窗 | |
| ''' | |
| def Barlett_window(length): | |
| res = [] | |
| for i in range(0, (length - 1) // 2, 1): # [0, (N-1)/2) | |
| res.append(2. * i / (length - 1)) | |
| for i in range((length - 1) // 2, N, 1): # [(N-1)/2, N-1] | |
| res.append(2. - 2. * i / (length - 1)) | |
| return res | |
| ''' | |
| 待滤波的sin(x) | |
| ''' | |
| def get_sin(dot_len): | |
| xs = np.linspace(-math.pi, math.pi, dot_len) | |
| ys = [math.sin(x) + random.uniform(-0.3, 0.3) for x in xs] | |
| return ys | |
| ''' | |
| 获取加窗结果 | |
| ''' | |
| def get_windowed(h_d, w_n): | |
| cutted = [] | |
| for i in range(N): | |
| cutted.append(h_d[i] * w_n[i]) | |
| return cutted | |
| ''' | |
| 把DFT结果X(k)插值到X(e^jw) | |
| ''' | |
| def interp(xk, dot_len=500): | |
| N = len(xk) | |
| res = [] | |
| eps = 0.000001 | |
| def phi(w): | |
| w += eps | |
| exp_part = Complex(math.cos((1 - N) * w / 2.), math.sin((1 - N) * w / 2.)) | |
| factor = Complex(1 / N * math.sin(N * w / 2.) / math.sin(w / 2.), 0) | |
| return factor * exp_part | |
| xs = np.linspace(0, math.pi, dot_len) | |
| for x in xs: | |
| sum = Complex(0, 0) | |
| for i in range(N): | |
| sum = sum + xk[i] * phi(x - 2. * math.pi * i / N) | |
| res.append(sum) | |
| return (xs, res) | |
| if __name__ == '__main__': | |
| # 生成滤波器无限长冲激响应并绘图 | |
| h_d_n = h_d(W_C) | |
| plt.title('滤波器的无限长冲激响应') | |
| plt.xlabel('n') | |
| plt.ylabel('h_d(n)') | |
| xs = range(0, N * 3, 1) | |
| plt.plot(xs, h_d_n, '.') | |
| plt.show() | |
| # 生成三种窗 | |
| R_N = rect_window(N) | |
| Barlett_N = Barlett_window(N) | |
| Hamming_N = Hamming_window(N) | |
| # 绘制三种加窗结果（时域） | |
| h_n_rect = get_windowed(h_d_n, R_N) | |
| h_n_Barlett = get_windowed(h_d_n, Barlett_N) | |
| h_n_Hamming = get_windowed(h_d_n, Hamming_N) | |
| plt.subplot(131) | |
| plt.xlabel('n') | |
| plt.ylabel('h(n)') | |
| plt.title('h(n) （矩形窗截断）') | |
| xs = range(0, N) | |
| plt.plot(xs, h_n_rect, '.') | |
| plt.subplot(132) | |
| plt.xlabel('n') | |
| plt.ylabel('h(n)') | |
| plt.title('h(n) （巴雷特窗截断）') | |
| plt.plot(xs, h_n_Barlett, '.') | |
| plt.subplot(133) | |
| plt.xlabel('n') | |
| plt.ylabel('h(n)') | |
| plt.title('h(n) （汉明窗截断）') | |
| plt.plot(xs, h_n_Hamming, '.') | |
| plt.show() | |
| # 将三种加窗转到频域，并绘制其幅度频谱 | |
| fft_h_n_rect = fft_dit2(convert_to_complex(add_zero_to_length(h_n_rect, 2 * N))) # 注意是2N点FFT | |
| fft_h_n_Barlett = fft_dit2(convert_to_complex(add_zero_to_length(h_n_Barlett, 2 * N))) | |
| fft_h_n_Hamming = fft_dit2(convert_to_complex(add_zero_to_length(h_n_Hamming, 2 * N))) | |
| plt.subplot(131) | |
| plt.xlabel('w') | |
| plt.ylabel('|H(e^jw)|') | |
| plt.title('|H(e^jw)|（矩形窗截断）') | |
| xs, dtft_h_n_rect = interp(fft_h_n_rect, 500) | |
| plt.plot(xs, convert_to_amplitude(dtft_h_n_rect), '-') | |
| plt.subplot(132) | |
| plt.xlabel('w') | |
| plt.ylabel('|H(e^jw)|') | |
| plt.title('|H(e^jw)|（巴雷特窗截断）') | |
| xs, dtft_h_n_Barlett = interp(fft_h_n_Barlett, 500) | |
| plt.plot(xs, convert_to_amplitude(dtft_h_n_Barlett), '-') | |
| plt.subplot(133) | |
| plt.xlabel('w') | |
| plt.ylabel('|H(e^jw)|') | |
| plt.title('|H(e^jw)|（汉明窗截断）') | |
| xs, dtft_h_n_Hamming = interp(fft_h_n_Hamming, 500) | |
| plt.plot(xs, convert_to_amplitude(dtft_h_n_Hamming), '-') | |
| plt.show() | |
| # 生成噪声干扰后的sin信号待滤波 | |
| sin_to_filter = get_sin(N) # 获取N点待滤波信号 | |
| plt.subplot(121) | |
| plt.xlabel('n') | |
| plt.ylabel('x(n)') | |
| plt.title('待滤波信号sin(n)+noise（时域）') | |
| plt.plot(range(0, N), sin_to_filter, '.') | |
| # 待滤波信号也做2N点FFT | |
| fft_sin_to_filter = fft_dit2(convert_to_complex(add_zero_to_length(sin_to_filter, 2 * N))) | |
| plt.subplot(122) | |
| plt.xlabel('k') | |
| plt.ylabel('X(k)') | |
| plt.title('待滤波信号sin(n)+noise（幅度频谱）') | |
| plt.plot(range(0, 2 * N), convert_to_amplitude(fft_sin_to_filter), '.') | |
| plt.show() | |
| # 使用汉明窗截断的滤波器进行低通滤波 | |
| fft_y_n_Hamming = [] | |
| for i in range(2 * N): | |
| fft_y_n_Hamming.append(fft_h_n_Hamming[i] * fft_sin_to_filter[i]) # 计算频域输出 | |
| plt.subplot(122) | |
| plt.xlabel('k') | |
| plt.ylabel('Y(k)') | |
| plt.title('滤波结果（幅度频谱）') | |
| plt.plot(range(0, 2 * N), convert_to_amplitude(fft_y_n_Hamming), '.') | |
| y_n_Hamming = convert_to_real(ifft(fft_y_n_Hamming)) # IFFT变回时域 | |
| y_n_Hamming = y_n_Hamming[0:-1] # 去掉最后一个，因为卷积结果长度应为2N-1 | |
| plt.subplot(121) | |
| plt.xlabel('n') | |
| plt.ylabel('y(n)') | |
| plt.title('滤波结果（时域）') | |
| plt.plot(range(0, 2 * N - 1, 1), y_n_Hamming, '.') | |
| plt.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment