一、FIR數字濾波器設計原理
本實驗采用窗函數法設計FIR數字低通濾波器。我們希望設計的濾波器系統函數如下：
\(H_{d}\left( e^{jw} \right) = \left\{ \begin{array}{l} {e^{- jw\alpha},~~~\left| w \right| \leq w_{c}} \\ {0,~~~{\rm otherwise}} \\ \end{array} \right.\)
它對應的單位沖激響應是：
\(h_{d}\left( n \right) = {\sin{\frac{ {w_{c} ( {n - \alpha} ) } }{\pi\left( n - \alpha \right)},~~~n \neq \alpha}}\)
\(w_c\)是截止頻率。
對它進行加窗操作后，單位沖激響應變為：
\(h\left( n \right) = h_{d}\left( n \right)w\left( n \right)\)
其中\(w(n)\)是窗函數的單位沖激響應。
為了滿足FIR濾波器的線性相位特性，我們取系統群時延：
\(\alpha = \frac {N-1}{2}\)
其中\(N\)是FIR濾波器的長度。這樣，\(h(n)\)可以關於\(n=α\)偶對稱。
二、基於FFT的濾波操作原理
記待濾波信號為\(x(n)\)，其長度也為\(N\)，則濾波結果信號：
\(y(n)=x(n)*h(n)\)
左右同作\(2N\)點FFT，有：
\({\rm FFT}_{2N} [y(n)]={\rm FFT}_{2N} [x(n)]∙{\rm FFT}_{2N} [h(n)]\)
作\(2N\)點FFT的原因是，\(y(n)\)的結果長度為\(2N-1\)，故至少應作\(2N-1\)點FFT才能獲得完整結果。又我們在實驗1實現的基2-FFT要求點數為2的冪次方，故應作\(2N\)點FFT，且實際上\(2N\)應是2的冪次。
再做反變換后有：
\(y(n)={\rm IFFT}[{\rm FFT}[y(n)]]\)
這樣就獲得了濾波結果信號的時域表示。
三、FIR數字濾波器的具體實現
程序首部可以調節\(w_c\)和\(N\)的取值。本報告中固定\(w_c=0.2π\)，\(N=32\)，后續不再說明。
首先設計濾波器的無限長沖激響應。我們利用\(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
接下來，我們用一個窗函數來截斷\(h_d (n)\)。窗函數的選取對濾波器的效果有很大影響。本實驗主要關注如下三種窗函數：
加窗相關代碼如下，以\(N\)為輸出長度作對應位相乘即可：
def get_windowed(h_d, w_n):
cutted = []
for i in range(N):
cutted.append(h_d[i] * w_n[i])
return cutted
以下是相關的時域圖像：
利用實驗1寫的相關代碼，並利用課本193頁（《數字信號處理教程》第四版，程佩青）從\(X(k)\)求取\(X(e^{jw})\)的插值公式，可以方便地繪制\(H(e^{jw})\)的幅度函數曲線，舉例代碼如下：
fft_h_n_rect = fft_dit2(convert_to_complex(add_zero_to_length(h_n_rect, 2 * N))) # 注意是2N點FFT
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), '-')
插值公式如下：
\(X(e^{{\rm j}w})=\sum_{k=0}^{N-1}{X(k)\Phi(\omega-2\pi k/N)}\)
\(\Phi(\omega)=\frac{1}{N} \frac{\sin(N\omega /2)}{\sin(\omega/2)}e^{-{\rm j}(N-1)\omega/2}\)
繪制結果如下：
從這三張幅度頻譜圖中已經可以清楚地看到不同的窗函數對濾波器的影響：
- 矩形窗的旁瓣抑制性能差，旁瓣峰值大，在阻帶仍有一定幅度值；巴特雷特窗旁瓣抑制性能稍好；漢明窗旁瓣抑制能力最強，在阻帶基本沒有幅度值。
- 矩形窗的主瓣窄；巴特雷特窗的主瓣寬度與漢明窗基本一致。
- 主瓣的寬度影響頻率的分辨率，旁瓣的強度影響干擾程度。在實際工程中，應在兩者之間進行權衡，選擇合適的窗函數。
四、使用FIR數字濾波器進行低通濾波
本節，我們利用上一節得到的漢明窗截斷的FIR濾波器進行低通濾波，使用被噪聲污染的正弦函數作為待濾波信號：
\(x\left( n \right) = {\sin n} + {\rm uniform\_ random}\left( - 0.3,~0.3 \right)\)
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
在\([-π,π]\)上采樣的\(N\)點\(x(n)\)的時域圖像和\(2N\)點幅度頻譜圖如下：
可以看到高頻分量具有較大的強度。
按照第二節提到的公式，在頻率域進行濾波，即對應位相乘：
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]) # 計算頻域輸出
然后把結果作反變換，獲得\(2N\)點時域序列，截去最后一個值，即為卷積的時域結果。
y_n_Hamming = convert_to_real(ifft(fft_y_n_Hamming)) # IFFT變回時域
y_n_Hamming = y_n_Hamming[0:-1] # 去掉最后一個，因為卷積結果長度應為2N-1
繪制相關圖像如下：
可以看到，高頻部分被有效抑制了，噪聲被有效消除了，正弦信號被較好地還原了。該FIR濾波器設計和應用成功。
附錄：完整代碼
fft.py請見實驗1
https://www.cnblogs.com/zxuuu/p/12425321.html
fir.py
# 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()