FIR連續(xù)采樣分段卷積時(shí)域重疊相加法

在上一個(gè)文檔里，我們提到了FIR系統(tǒng)在時(shí)域的分段卷積中使用“重疊保留(Overlap-Save)”的處理方式，這里我們續(xù)集，說(shuō)明一下“重疊相加(Overlap-Add)”的處理方式。 信號(hào)處理在時(shí)域和頻域中處理是有差異的。 說(shuō)得通俗一點(diǎn)就是：時(shí)域中處理是直接用采集到的信號(hào)進(jìn)行計(jì)算；而頻域中則要用離散傅里葉變換(DFT)/離散傅里葉反變換(IDFT)對(duì)采集到的信號(hào)進(jìn)行轉(zhuǎn)換到頻域，然后再?gòu)念l域轉(zhuǎn)回時(shí)域處理。

在我們看到的參考文檔[1]中，描述的是在頻域進(jìn)行DFT，然后由IDFT轉(zhuǎn)回時(shí)域處理的過(guò)程。如下圖所示。大概的過(guò)程是：

每次處理的數(shù)據(jù)長(zhǎng)度為L(zhǎng)，然后在每分段的尾部添加（M-1）個(gè)0之后，讓每次處理的數(shù)據(jù)序列長(zhǎng)度N=L+M-1，通常N為2的冪次倍；

同時(shí)對(duì)于濾波器，也需要將原來(lái)長(zhǎng)度為b的序列，通過(guò)填0的方式增加到長(zhǎng)度為N；

由DFT將兩個(gè)數(shù)列分別轉(zhuǎn)換到頻域，相乘后，再I(mǎi)DFT轉(zhuǎn)回時(shí)域，就得到N=（L+M-1）的時(shí)域卷積結(jié)果；

保留每次操作的所有數(shù)據(jù)，然后在下一次操作結(jié)束后，將最新數(shù)據(jù)的最前面的（M-1）個(gè)結(jié)果數(shù)據(jù)和上一次結(jié)果數(shù)據(jù)的最后（M-1）個(gè)數(shù)據(jù)順序相加......持續(xù)直至結(jié)束。

FIR頻域的重疊相加示意圖

我們看看時(shí)域的卷積應(yīng)該怎么操作。

FIR時(shí)域重疊相加操作示意圖

如上圖所示：

每次讀取長(zhǎng)為L(zhǎng)的數(shù)據(jù)序列，然后與長(zhǎng)度為M的濾波器進(jìn)行卷積，生成一個(gè)（L+M-1）的卷積結(jié)果序列；

保留每次操作的所有數(shù)據(jù)，然后在下一次操作結(jié)束后，將最新數(shù)據(jù)的最前面的（M-1）個(gè)結(jié)果數(shù)據(jù)和上一次結(jié)果數(shù)據(jù)的最后（M-1）個(gè)數(shù)據(jù)順序相加......持續(xù)直至結(jié)束。

兩個(gè)過(guò)程看起來(lái)略有差異，甚至?xí)X(jué)得時(shí)域的處理更簡(jiǎn)單省事，會(huì)不會(huì)更省時(shí)？其實(shí)頻域看起來(lái)費(fèi)時(shí)，但數(shù)據(jù)規(guī)模到了一定程度之后，頻域的處理速度就具有優(yōu)勢(shì)了。然而對(duì)于一般的應(yīng)用，直接卷積操作還是可以接受的。

還記得上次我們提到Python中卷積函數(shù)np.convolve的三種模式吧？該函數(shù)對(duì)卷積是在時(shí)域中進(jìn)行的。

在Python中，卷積函數(shù)np.convolve(data_segment, b, mode)對(duì)指定長(zhǎng)度的數(shù)據(jù)data_segment（長(zhǎng)度L），和FIR濾波器系數(shù)序列b（長(zhǎng)度M）進(jìn)行卷積。輸出的結(jié)果序列則分為以下三種：

full：結(jié)果長(zhǎng)度=M+L-1

same:結(jié)果長(zhǎng)度=max(M,L)

valid:結(jié)果長(zhǎng)度=max(M,L)-min(M,L)+1=L-(M-1)

在這里，我們需要選用full模式，這樣就獲取每段卷積一個(gè)不落的所有數(shù)據(jù)（L+M-1)。先看模擬效果后看Python代碼。 故事情節(jié)設(shè)定：50Hz的信號(hào)中，夾雜300，450Hz的干擾。濾除干擾。

FIR選頻濾波器的幅頻響應(yīng)

FIR系統(tǒng)重疊相加的濾波結(jié)果示意圖

這里要特別說(shuō)明一點(diǎn)：卷積后的數(shù)據(jù)長(zhǎng)度，最終會(huì)比原來(lái)的數(shù)多出（M-1）個(gè)，所以輸出到圖的時(shí)候，需要有意控制長(zhǎng)度。 濾波過(guò)程中要經(jīng)歷“熱身”，所以最開(kāi)始階段有（M-1）個(gè)數(shù)據(jù)也是可以剔除的。同樣，如果我們看卷積最終結(jié)果尾部不處理，也有（M-1）個(gè)無(wú)效數(shù)據(jù)的輸出需要截取。

卷積后尾部無(wú)效數(shù)據(jù)（M-1）

上代碼，我們自己在代碼中劃重點(diǎn)，并調(diào)整輸出結(jié)果的有效長(zhǎng)度范圍：

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fft import fft
import math


# 創(chuàng)建帶通濾波器
f1 = 40
f2 = 60
filter_len = 80     # 濾波器長(zhǎng)度
fs = 1600           # 采樣頻率維持不變


b = signal.firwin(filter_len, [f1, f2], pass_zero=False, fs=fs)


# 設(shè)置數(shù)據(jù)長(zhǎng)度
seg_filter_len = 256                            # filter output length of each segment data
segment_len = seg_filter_len - filter_len + 1   # 分段數(shù)據(jù)目標(biāo)長(zhǎng)度 seg_filter_len = segment_len + filter_len - 1
target_length = segment_len * 50                # 總數(shù)據(jù)長(zhǎng)度


# 而新的時(shí)間序列的上限 b
bspace = target_length / fs


# 生成的時(shí)間序列為L(zhǎng)的整數(shù)倍，模擬每次采樣的數(shù)據(jù)的長(zhǎng)度
t = np.arange(0, bspace, 1/fs)


# 產(chǎn)生一個(gè)含有300Hz,450Hz和50Hz信號(hào)的模擬信號(hào)
x = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 300 * t) + 0.5 * np.sin(2 * np.pi * 450 * t)


segments = []
for i in range(0, len(x), segment_len):
    segments.append(x[i:i+segment_len])


#Filtering&Overlap-Add processing
# Total outputput buffer, len = target_length + filter_len - 1
filtered_signals = np.zeros(target_length + filter_len - 1)             
for i in range(len(segments)):
    filtered_segment = np.convolve(segments[i], b, mode='full')    # full模式用于保留所有卷積結(jié)果 N = L + M -1
    filtered_signals[i*segment_len:i*segment_len+len(filtered_segment)] += filtered_segment # 疊加過(guò)程


filtered_signals = filtered_signals#[:target_length] # 保留和原信號(hào)等長(zhǎng)的部分


#FilterFreqResponse
w, h = signal.freqz(b, 1, fs=fs)
plt.figure()
plt.plot(w, abs(h))
plt.title('Filter Freq Response')
plt.grid()
plt.xlabel('f[Hz]')
plt.ylabel('Amplitude')


# Signal Before filtering & Spectrum
n = len(x)
freq = np.fft.fftfreq(n, 1/fs)
y = np.fft.fft(x)


plt.figure()
plt.subplot(221)
plt.plot(t[:500], x[:500])
plt.title('Original Signal')
plt.subplot(222)
plt.plot(freq[:n//2], np.abs(y[:n//2]*2/n)) # 標(biāo)幺，繪制前一半
plt.title('SpectrumofOrginalSignal')
plt.grid()


#SignalAfterfiltering & Spectrum
n = len(x)
y = np.fft.fft(filtered_signals)


plt.subplot(223)
# 1. Normal output
plt.plot(t[:500], filtered_signals[:500])
plt.title('Filtered Signal')
plt.subplot(224)
plt.plot(freq[:n//2], np.abs(y[:n//2]*2/n)) # 標(biāo)幺，繪制前一半
plt.title('Spectrum of Filtered Signal')
plt.grid()


plt.tight_layout()


plt.show()


# 2. End of convolution without end cutoff: for test purpose
plt.figure()
temp = t[8000:8850]
s = filtered_signals[8079:8929]
plt.plot(temp, s)
plt.title('Filtered Signal')
plt.grid()
plt.show()