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

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

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

同時對于濾波器，也需要將原來長度為b的序列，通過填0的方式增加到長度為N；

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

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

FIR頻域的重疊相加示意圖

我們看看時域的卷積應該怎么操作。

FIR時域重疊相加操作示意圖

如上圖所示：

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

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

兩個過程看起來略有差異，甚至會覺得時域的處理更簡單省事，會不會更省時？其實頻域看起來費時，但數(shù)據(jù)規(guī)模到了一定程度之后，頻域的處理速度就具有優(yōu)勢了。然而對于一般的應用，直接卷積操作還是可以接受的。

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

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

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

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

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

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

FIR選頻濾波器的幅頻響應

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

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

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

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

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     # 濾波器長度
fs = 1600           # 采樣頻率維持不變


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


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


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


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


# 產(chǎn)生一個含有300Hz,450Hz和50Hz信號的模擬信號
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 # 疊加過程


filtered_signals = filtered_signals#[:target_length] # 保留和原信號等長的部分


#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)) # 標幺，繪制前一半
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)) # 標幺，繪制前一半
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()