在之前我們的博文中，所提到的數(shù)據(jù)的濾波處理和仿真分析，其實(shí)都是圍著一段固定長(zhǎng)度的模擬數(shù)據(jù)展開(kāi)的，除了知道濾波器的幅頻、相頻響應(yīng)特性之外，也直觀地看到了濾波的效果會(huì)是怎么樣的。實(shí)際應(yīng)用會(huì)怎么樣？需要怎么處理？

假設(shè)需要對(duì)音頻，視頻，或者儲(chǔ)存的大量數(shù)據(jù)等進(jìn)行處理時(shí)，期待全部讀取然后一次處理顯然是不合適，也沒(méi)有必要，而且資源和實(shí)時(shí)性都無(wú)法滿足。而分段處理，需要的資源少是一方面，并行和分布處理難道不考慮一下？

那分段處理，能否將每段濾波之后獲取的數(shù)據(jù)直接拼接作為整體濾波后的數(shù)據(jù)呢？然而分段處理需要考慮有其特殊性。

就信號(hào)處理而言，F(xiàn)IR系統(tǒng)和IIR系統(tǒng)的處理方式各有特點(diǎn)。

IIR系統(tǒng)的連續(xù)采集濾波

IIR系統(tǒng)的輸出不僅依賴于當(dāng)前和過(guò)去的激勵(lì)輸入，而且還依賴于系統(tǒng)過(guò)去的輸出（反饋）。IIR系統(tǒng)的特性是有無(wú)限長(zhǎng)的沖擊響應(yīng)，理論上即使停止了前饋輸入，它還有系統(tǒng)過(guò)去的輸出繼續(xù)作為反饋輸入，從而維持輸出。如果一定要卷積IIR系統(tǒng)，通常要做適當(dāng)?shù)慕財(cái)嗪徒铺幚恚趴梢宰孖IR濾波的效果在有限的計(jì)算資源使用卷積。不過(guò)雖然理論可卷，實(shí)際多用時(shí)域遞推的方式。

IIR系統(tǒng)的遞推計(jì)算，之前我們提供的示例C++濾波器代碼中，是為了演示，一次性批發(fā)處理多個(gè)數(shù)據(jù)，實(shí)際處理是可以根據(jù)應(yīng)用來(lái)調(diào)整每次輸入數(shù)據(jù)的長(zhǎng)度。無(wú)論是分段處理還是單個(gè)輸出，如果保留了過(guò)去的輸入和輸出狀態(tài)，是可以持續(xù)進(jìn)行處理的。下面是每次輸入一個(gè)信號(hào)x(n)，然后生成/讀取一個(gè)輸出y(n)，一定要正確保留每次的輸出和輸入。

下面是一個(gè)簡(jiǎn)單IIR單輸入濾波的示例，參數(shù)沒(méi)有哈。當(dāng)然，把這個(gè)轉(zhuǎn)換為FIR的濾波處理也是很自然的事情。

#define N 3          // 設(shè)置濾波器的階數(shù)


double a[N] = {...}; // 這是濾波器的反饋系數(shù)
double b[N] = {...}; // 這是濾波器的前饋系數(shù)
double x[N] = {0};   // 這是輸入的緩沖數(shù)組
double y[N] = {0};   // 這是輸出的緩沖數(shù)組


double iir_filter(double input)
{
    int i;
//更新輸入緩沖數(shù)組，循環(huán)更替
    for (i = N - 1; i > 0; --i)
    {
        x[i] = x[i - 1];
    }
    x[0] = input;
//計(jì)算當(dāng)前的輸出
    double output = 0;
    for (i = 0; i < N; ++i)
    {
        output += b[i] * x[i] - a[i] * y[i];
    }
????//?更新輸出緩沖數(shù)組
    for (i = N - 1; i > 0; --i)
    {
        y[i] = y[i - 1];
    }
    y[0] = output;
returnoutput;
}

如果IIR系統(tǒng)長(zhǎng)數(shù)據(jù)分段沒(méi)有保存過(guò)去的狀態(tài)，我們可以得到下圖。

IIR分段濾波沒(méi)有保留過(guò)去的狀態(tài)

過(guò)濾沒(méi)有處理好，還生成了額外的干擾信號(hào)。總之，IIR系統(tǒng)中忘記歷史，就會(huì)有懲戒啊。

如果IIR系統(tǒng)長(zhǎng)數(shù)據(jù)分段保存過(guò)去的狀態(tài)，我們可以得到下圖：

IIR分段濾波有過(guò)去的狀態(tài)

前面兩組圖相比，差異是不是一目了然？下面是IIR系統(tǒng)每次分段考慮了過(guò)去狀態(tài)的Python模擬代碼。

# IIR 系統(tǒng)分段處理：過(guò)去狀態(tài)的處理
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from numpy.fft import fft, fftfreq


# 創(chuàng)建模擬信號(hào)
# 采樣頻率
Fs = 1000
nyq = 0.5 * Fs


# 設(shè)定阻帶頻率
low = 45.0
high = 55.0


# 按照通帶和阻帶頻率的順序來(lái)設(shè)計(jì)濾波器
lowcut = low / nyq
highcut = high / nyq


# 使用iirdesign設(shè)計(jì)濾波器
b, a = signal.iirdesign([lowcut, highcut], [lowcut - 0.05, highcut + 0.05], gpass=1, gstop=60, ftype='butter', output='ba')


T = 1.0 / Fs # 采樣時(shí)間間隔
t = np.arange(0, 1, T) # 創(chuàng)建時(shí)間軸
x = np.sin(2 * np.pi * 50 * t) + 1 * np.sin(2 * np.pi * 220 * t) # 50Hz和120Hz信號(hào)疊加
# 創(chuàng)建一個(gè)與x長(zhǎng)度相同的用于保存濾波后結(jié)果的數(shù)組
y = np.zeros_like(x)


# 分段進(jìn)行濾波
win_size = 100 # 窗口大小
zi = signal.lfilter_zi(b, a) * x[0] # 計(jì)算初始狀態(tài)
for i in range(win_size, len(x)+1, win_size):
    y[i-win_size:i], zi = signal.lfilter(b, a, x[i-win_size:i], zi=zi)


# 計(jì)算輸入信號(hào)和輸出信號(hào)的頻譜
frequencies = fftfreq(len(t), 1/Fs)
x_spectrum = np.abs(fft(x))
y_spectrum = np.abs(fft(y))


# 繪制幅頻響應(yīng)圖
w, h = signal.freqz(b, a)
plt.figure()
plt.plot(w/(2*np.pi), abs(h))
plt.title('Frequency response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Gain')
plt.grid()


# 繪制頻譜圖
plt.figure(figsize=(9, 6))
plt.subplot(221)
plt.title('Spectrum Before filtering')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')
plt.grid()
plt.plot(frequencies[:Fs//2], x_spectrum[:Fs//2], label='Before filtering')
plt.subplot(222)
plt.plot(frequencies[:Fs//2], y_spectrum[:Fs//2], label='After filtering')
plt.title('Spectrum After filtering')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')
plt.grid()


# 繪制模擬信號(hào)和濾波后的信號(hào)
plt.subplot(223)
plt.plot(t, x, label='Before filtering')
plt.title('Time domain Before filtering')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.grid()
plt.subplot(224)
plt.plot(t, y, label='After filtering')
plt.title('Time domain After filtering')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.grid()
plt.tight_layout()
plt.legend()


plt.show()

FIR系統(tǒng)的連續(xù)采集濾波

和IIR系統(tǒng)相比，F(xiàn)IR系統(tǒng)的輸出則僅僅依靠當(dāng)前和過(guò)去的輸入。但是，要實(shí)現(xiàn)相同的濾波效果（不考慮相位線性），F(xiàn)IR往往要更高的階數(shù)才可以完成。階數(shù)越高，意味著越多的資源和越長(zhǎng)的系統(tǒng)延遲。

我們知道，F(xiàn)IR濾波器中的每個(gè)系數(shù)b(n)，都對(duì)應(yīng)于系統(tǒng)的單位沖擊響應(yīng)h(n)。所以，對(duì)于FIR濾波器，那我們就可以通過(guò)但不限于卷積的方式來(lái)處理信號(hào)。這里我們圍繞卷積的方式來(lái)討論FIR濾波器處理數(shù)據(jù)分段濾波的兩種方法：

重疊保留(Overlap-Save)

重疊相加(Overlap-Add)

下面兩個(gè)等式分別是卷積和FIR濾波處理。

（1）

（2）

如果把b[k]看作h(k)，則水到渠成，F(xiàn)IR濾波器的處理過(guò)程就是一個(gè)卷積運(yùn)算。變成卷積后的另外的好處就是，我們可以利用FFT/IFFT來(lái)進(jìn)行某些大型數(shù)據(jù)的加速運(yùn)算——時(shí)域的卷積等同于頻域的乘積。對(duì)于小到中等大小的數(shù)組序列，一般也就用時(shí)域的卷積了。

但是卷積從現(xiàn)象上看，當(dāng)我們參與卷積的FIR濾波器長(zhǎng)度為M，被卷積的數(shù)據(jù)序列長(zhǎng)度為L(zhǎng)，那么，每次卷積的結(jié)果的數(shù)據(jù)序列N長(zhǎng)度會(huì)是：

N= M + L -1

當(dāng)你面對(duì)每次多出原先數(shù)據(jù)序列長(zhǎng)度L的結(jié)果N時(shí)？如何取舍？

如果濾波器是從全0的狀態(tài)開(kāi)始的，濾波器需要一定的"預(yù)熱"時(shí)間才能達(dá)到穩(wěn)定的工作狀態(tài)。濾波器的"預(yù)熱"過(guò)程就是濾波器讀取前N個(gè)樣本的過(guò)程。因此，你會(huì)看到濾波后的數(shù)據(jù)最開(kāi)始部分是從0開(kāi)始，然后逐漸看到正常的輸出波形。

另外，我們示例的濾波器是加窗處理的，而且模擬信號(hào)頻率比較單一，相距也較遠(yuǎn)，否則，處理不好的濾波結(jié)果會(huì)讓人感覺(jué)出乎意外。因?yàn)槊看谓厝∫欢螖?shù)據(jù)，相對(duì)于給采集的數(shù)據(jù)加了一個(gè)矩形窗，而我們?cè)谠O(shè)計(jì)濾波器時(shí)，是需要對(duì)這個(gè)矩形窗進(jìn)行額外加窗處理的。矩形窗對(duì)于波形而言存在著頻譜泄漏和主瓣寬度等問(wèn)題，對(duì)數(shù)據(jù)的處理有時(shí)候會(huì)比較麻煩。——注意，要跑題了吧？

我們看模擬結(jié)果。下圖為待濾波處理的模擬信號(hào)。500Hz+1400Hz。

下圖為矩形窗濾波器的幅頻響應(yīng)——知道為什么我要選1400Hz作為被濾掉的信號(hào)了吧？

矩形窗選頻濾波器幅頻特性

這個(gè)例子僅僅是為了說(shuō)明矩形窗的影響。下面是矩形窗濾波器的結(jié)果，500Hz+1400Hz的信號(hào)經(jīng)過(guò)濾波之后的輸出。——跑題十萬(wàn)八千里是為了一扇窗？

（忽略此圖）

矩形窗濾波器處理結(jié)果

而加漢明窗的濾波結(jié)果在試驗(yàn)設(shè)定條件下則要好一些，同樣留意1400Hz處。如下圖。

漢明窗濾波器的幅頻特性

加漢明窗的濾波結(jié)果

以上處理為了差異效果，對(duì)于濾波器長(zhǎng)度、一次處理數(shù)據(jù)長(zhǎng)度，兩個(gè)信號(hào)頻率的距離都是作了特別的處理。誰(shuí)說(shuō)實(shí)際應(yīng)用中不會(huì)碰到呢？

帶通濾波器不同的加窗幅頻特性

以下內(nèi)容，都默認(rèn)已經(jīng)加窗。

我們?cè)俅位氐紽IR系統(tǒng)使用重疊保留(Overlap-Save)和重疊相加(Overlap-Add)的方式進(jìn)行分段處理這個(gè)話題。

如果FIR系統(tǒng)在頻域通過(guò)DFT/IDFT并采用重疊保留的方式，那應(yīng)該是按照下圖的方式進(jìn)行的。

對(duì)于FIR濾波器長(zhǎng)度為M，數(shù)據(jù)段長(zhǎng)度為L(zhǎng)時(shí)，流程大致是：在第一個(gè)數(shù)據(jù)段加M-1個(gè)0，其余的數(shù)據(jù)段都是加上前一段的最后M-1個(gè)數(shù)據(jù)；然后進(jìn)行DFT/IDFT處理得到長(zhǎng)度為(L+M-1)的y(n)序列，然后去掉前面M-1個(gè)無(wú)效值，保留后面的L個(gè)。

頻域中FIR系統(tǒng)的重疊保留法示意圖

在時(shí)域我們也想照貓畫(huà)虎的。然而，對(duì)于FIR濾波器長(zhǎng)度為M和信號(hào)的長(zhǎng)度L進(jìn)行卷積，其結(jié)果的長(zhǎng)度將會(huì)是M + L - 1。在實(shí)際應(yīng)用中，我們通常關(guān)注的部分是兩個(gè)序列完全重疊的區(qū)域，也就是“有效值”。這意味著如果兩個(gè)參與卷積的長(zhǎng)度不一樣（一般L>M)，那最后的“有效值”是需要在卷積的結(jié)果序列上兩頭都要掐掉（M-1）個(gè)元素。所以實(shí)際卷積結(jié)果中的有效序列長(zhǎng)度是：

[L+M-1-2（M-1)]=L-M+1

在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)

為了測(cè)試時(shí)域中的重疊保留處理方式，小編選了上面卷積函數(shù)的valid模式。

在重疊保留方式下，每次參與卷積的數(shù)據(jù)長(zhǎng)度不是L，而是每個(gè)數(shù)據(jù)L起始處，都額外添加了（M-1）長(zhǎng)的0（第一個(gè)）或者前一段數(shù)據(jù)的最后（M-1）個(gè)數(shù)據(jù)，再與長(zhǎng)度為M的濾波器卷積。這種情況下，參與卷積的數(shù)據(jù)長(zhǎng)度和生成的卷積結(jié)果長(zhǎng)度是有變化的。

L和M長(zhǎng)序列進(jìn)行卷積得到的full長(zhǎng)度為：M+L-1；

那如果將L+(M-1)和M長(zhǎng)的序列進(jìn)行卷積，得到的full長(zhǎng)度為：L+(M-1)+M-1=L+2(M-1)；

根據(jù)步驟2，我們把得到的full長(zhǎng)度的結(jié)果兩邊減去(M-1)就可以得到一個(gè)長(zhǎng)度為L的有效卷積序列值。

所以這個(gè)操作和頻域中DFT、IDFT之后得到的保留方式有差異，因?yàn)闀r(shí)域中卷積之后的分段結(jié)果數(shù)據(jù)的長(zhǎng)度會(huì)由L+M-1變?yōu)長(zhǎng)+2(M-1)，最終取卷積有效值段的長(zhǎng)度就等于[L+2(M-1)-2(M-1)]=L。

下面提供模擬仿真。輸出圖和Python模擬代碼。（沒(méi)有畫(huà)時(shí)域卷積中FIR系統(tǒng)的數(shù)據(jù)長(zhǎng)度示意圖，無(wú)他，偷懶。）

在500Hz信號(hào)中混有1400Hz和2400Hz的干擾。采用分段濾波的方式對(duì)FIR濾波器進(jìn)行卷積處理。我們?nèi)匀话阉麨椤爸丿B保留法”。

FIR系統(tǒng)時(shí)域卷積的重疊保留法

# FIR系統(tǒng)時(shí)域卷積：重疊保留法
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fft import fft,fftfreq


# 設(shè)置采樣率和數(shù)據(jù)長(zhǎng)度
fs = 6000


# 創(chuàng)建帶通濾波器
# Create a bandpass filter
f1 = 400
f2 = 600
filter_len = 100         # FIR濾波器長(zhǎng)度
b = signal.firwin(filter_len, [f1, f2], pass_zero=False, fs=fs, window='hamming')  #boxcar 


# 設(shè)置數(shù)據(jù)長(zhǎng)度：模擬每次采集同樣長(zhǎng)度的數(shù)據(jù)
# segment_len，是按照np.concatenate函數(shù)在Valid模式下卷積后有效長(zhǎng)度來(lái)計(jì)算的
seg_filter_len = 512                            # 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í)間序列的上限
TimeSpace = target_length / fs


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


#產(chǎn)生一個(gè)含有1.4KHz，2.4KHz和500Hz信號(hào)的模擬信號(hào)，其中1.4，2.4KHz的信號(hào)將被過(guò)濾
x = np.sin(2 * np.pi * 500 * t) + 0.5 * np.sin(2 * np.pi * 1400 * t) + 0.5 * np.sin(2 * np.pi * 2400 * t)


# Filter and segment with overlap-save method
# each segment length should be longer than filter's length
y_result = []
segment = None
for i in range(0, target_length, segment_len):
    if i==0:
        # 第一段數(shù)據(jù)濾波前，前面填充(M-1)個(gè)0，再加開(kāi)始的L個(gè)數(shù)據(jù)參與濾波卷積
        segment = np.concatenate((np.zeros(filter_len-1), x[0:segment_len])) 
    else:
        # 對(duì)后續(xù)的每個(gè)段，包括從上段的末尾取(M-1)個(gè)點(diǎn)，再新取L個(gè)點(diǎn)，這里M指濾波器長(zhǎng)度,L指擬當(dāng)前新取的數(shù)據(jù)長(zhǎng)度
        segment = np.concatenate((x[i-filter_len+1:i], x[i:i+segment_len]))
    
    # 每個(gè)段使用“valid”卷積模式，所以每次卷積后的有效數(shù)據(jù)長(zhǎng)度是:L = L + 2(M-1) - 2(M-1) = segment_len
    conv = np.convolve(segment, b, mode='valid')
    # 與結(jié)果進(jìn)行合并連接
    y_result = np.concatenate((y_result, conv))


# Frequency response of filter
w, h = signal.freqz(b, 1, fs=fs)
plt.figure()
plt.plot(w, abs(h))
plt.title('Frequency Response')     # 頻率響應(yīng)
plt.xlabel('Frequency [Hz]')        # 頻率
plt.ylabel('Amplitude')             # 幅度
plt.grid()


# Plot the original signal and its FFT
n = len(x)
freq = fftfreq(n, 1/fs)
y = fft(x)


plt.figure(figsize=(9,6))
plt.subplot(221)
plt.plot(t[:500], x[:500])
plt.title('Original Signal')                # 原始信號(hào)
plt.grid()
plt.subplot(222)
plt.plot(freq[:n//2], np.abs(y[:n//2]*2/n)) # 頻譜規(guī)范化輸出
plt.title('FFT of Original Signal')         # 原始信號(hào)的FFT
plt.grid()


# Plot the filtered signal and its FFT
n = len(y_result)
freq = fftfreq(n, 1/fs)
y = fft(y_result)


plt.subplot(223)
plt.plot(t[:500], y_result[:500])
plt.title('Filtered Signal')                # 濾波后的信號(hào)
plt.grid()
plt.subplot(224)
plt.plot(freq[:n//2], np.abs(y[:n//2]*2/n)) # 頻譜規(guī)范化輸出
plt.title('FFT of Filtered Signal')         # 濾波后信號(hào)的FFT
plt.grid()


plt.tight_layout()
plt.show()

大膽假設(shè)，小心求證。

以上代碼并未實(shí)證，純屬模擬，各位有興趣可以試試。

不過(guò)我們安費(fèi)諾的傳感器都是久經(jīng)沙場(chǎng)，經(jīng)歷過(guò)各種工況驗(yàn)證的，感興趣也可以試試。

內(nèi)容長(zhǎng)了看起來(lái)都煩，F(xiàn)IR時(shí)域的“重疊相加”卷積處理，下回分解。