# 自適應濾波器和NLMS模擬 本章將簡要介紹自適應濾波器的原理以及其最常用的算法NLMS，并給NLMS算法的兩種實現方法：用純Python編寫，和用ctypes調用C語言編寫。最后將對NLMS算法進行一些的實驗。 ## 自適應濾波器簡介 近年來，隨著數字信號處理器的功能的不斷增強，自適應信號處理 (adaptive signal process)活躍在噪聲消除、回聲控制、信號預測、聲音定位等眾多信號處理領域。 盡管其應用領域十分廣泛，但基本的系統構造大致只有如下幾種分類。 ### 系統辨識 所謂系統辨識(system identification)，就是通過對未知系統的輸入輸出進行觀測，構造一個濾波器使得它在同樣的輸入的情況下，輸出信號和未知系統相同。簡而言之，就是通過觀測未知系統對輸入的反應，探知其內部情況。為了探知內情而使用的輸入信號我們稱之為參照信號。  系統識別(System Identification)的框圖 如上圖所示參照信號 x(j)同時輸入到未知系統和自適應濾波器H中，未知系統的輸出為y(j)， 自適應濾波器的輸出為u(j)，由于觀測誤差或者外部噪聲的干擾，實際觀測到的未知系統的輸出為d(j)=y(j)+n(j)，n(j)被稱為外部干擾。通過求的d(j)和u(j)之間的誤差e(j)=d(j)-u(j)，我們可以知道自適應濾波器H和未知系統還有多少差別，通過這個誤差我們更新H的內部參數，使得它更加靠近未知系統。 上面各個公式中的j表示某一時刻，因為我們討論的是數字信號處理，已經對所有的信號進行取樣，因此可以把j簡單的看作取樣點的下標。 ### 信號預測 所謂信號預測就是通過信號過去的值預測（計算）現在的值，下面是信號預測的系統框圖。  信號預測(Predication)框圖 x(j)是待預測的信號，假設我們無法完美地觀測此信號，因此導入一個外部干擾n(j)，這樣d(j)=x(j)+n(j)就是我們觀測到的待預測信號。 通過延遲器將d(j)進行延時得到d(j-D)，并把d(j-D)輸入到自適應濾波器H中，得到其輸出為u(j)，u(j)就是自適應濾波器通過待預測信號過去的值預測出的現在的值，計算觀測值d(j)和預測值u(j)之間的誤差e(j)=d(j)-u(j)，通過e(j)更新自適應濾波器H的內部系數使得其輸出更加接近d(j)。 如果x(j)存在白色噪聲的成分和周期信號的成分，由于白色噪聲是完全不自相關，無法預測的信號，因此通過過去的值x(j-D)所能預測的只能是其中的周期信號的成分。這樣自適應濾波器H的輸出信號u(j)就會與周期信號成分漸漸逼近，而e(j)則是剩下的不可預測的白色噪聲的成分。因此自適應濾波器也可以運用于噪聲消除。 ### 信號均衡  信號均衡(Equalization)框圖 當信號x(j)通過未知系統之后變成y(j)，未知系統對信號x(j)進行了某種改變，使得其波形產生歪曲。我們希望均衡器矯正這種歪曲，也就是通過y(j)重建原始信號x(j)，由于因果律還原原始信號x(j)是不可能的，我們只能還原其延時了的信號x(j-D)。x(j)和x(j-D)除了時間上的延遲之外，其它特性完全相同。 這里我們將觀測到的未知系統的輸出y(j)+n(j)輸入到自適應濾波器H中，通過H的系數更新使得其輸出u(j)逐漸逼近原始信號的延時x(j-D)。這樣我們就構建了一個濾波器H使得它與未知系統的卷積正好等于一個脈沖傳遞函數。也就是說H的頻域特性恰好能抵消未知系統的所帶來的改變。 ## NLMS計算公式 自適應濾波器中最重要的一個環節就是其系數的更新算法，如果不對自適應濾波器的系數更新的話，那么它就只是一個普通的FIR濾波器了。系數更新算法有很多種類，最基本、常用、簡單的一種方法叫做NLMS(歸一化最小均方)，讓我們先來看看它的數學公式表達： 設置自適應濾波器系數  的所有初始值為0,  的長度為I。  對每個取樣值進行如下計算，其中n=0, 1, 2, ...    自適應濾波器系數  是一個長度為I的矢量，也就是一個長度為I的FIR濾波器。在時刻n，濾波器的每個系數對應的輸入信號為  ，它也是一個長度為I的矢量。這兩個矢量的點乘即為濾波器的輸出和目標信號d(n)之間的差為e(n)，然后根據e(n)和  ， 更新濾波器的系數。 數學公式總是令人難以理解的，下面我們以圖示為例進行說明：  NLMS算法示意圖 圖中假設自適應濾波器h的長度為4，在時刻7濾波器的輸出為: ``` u[7] = h[0]*x[7] + h[1]*x[6] + h[2]*x[5] + h[3]*x[4] ``` 濾波器的輸入信號的平方和powerX為: ``` powerX = x[4]*x[4] + x[5]*x[5] + x[6]*x[6] + x[7]*x[7] ``` 未知系統的輸出d[7]和濾波器的輸出u[7]之間的差為: ``` e[7] = d[7] - u[7] ``` 使用u[7]和x[4]..x[7]對濾波器的系數更新: ``` h[4] = h[4] + u * e[7]*x[4]/powerX h[4] = h[5] + u * e[7]*x[5]/powerX h[4] = h[6] + u * e[7]*x[6]/powerX h[4] = h[7] + u * e[7]*x[7]/powerX ``` 其中參數u成為更新系數，為0到1之間的一個實數，此值越大系數更新的速度越快。對于每個時刻i都需要進行上述的計算，因此濾波器的系數對于每個參照信號x的取樣都更新一次。 ## NumPy實現 按照上面介紹的NLMS算法，我們很容易寫出用NumPy實現的NLMS計算程序： ``` # -*- coding: utf-8 -*- # filename: nlms_numpy.py import numpy as np # 用Numpy實現的NLMS算法 # x為參照信號，d為目標信號，h為自適應濾波器的初值 # step_size為更新系數 def nlms(x, d, h, step_size=0.5): i = len(h) size = len(x) # 計算輸入到h中的參照信號的乘方he power = np.sum( x[i:i-len(h):-1] * x[i:i-len(h):-1] ) u = np.zeros(size, dtype=np.float64) while True: x_input = x[i:i-len(h):-1] u[i] = np.dot(x_input , h) e = d[i] - u[i] h += step_size * e / power * x_input power -= x_input[-1] * x_input[-1] # 減去最早的取樣 i+=1 if i >= size: return u power += x[i] * x[i] # 增加最新的取樣 ``` 為了節省計算時間，我們用一個臨時變量power保存輸入到濾波器h中的參照信號x的能量。在對于x中的每個取樣的循環中，power減去x中最早的一個取樣值的乘方，增加最新的一個取樣值的乘方。這樣為了計算參照信號的能量，每次循環只需要計算兩次乘法和兩次加法即可。 nlms函數的輸入為參照信號x、目標信號d和自適應濾波器的系數h。因為在后面的模擬計算中，d是x和未知系統的脈沖響應的卷積而計算的來，它的長度會大于x的參數，因此循環體的循環次數以參照信號的長度為基準。 為了對自適應濾波器的各種應用進行模擬，我們還需要如下的幾個輔助函數，完整的程序請參考 [_NLMS算法的模擬測試_](example_nlms_test.html) 。 ``` def make_path(delay, length): path_length = length - delay h = np.zeros(length, np.float64) h[delay:] = np.random.standard_normal(path_length) * np.exp( np.linspace(0, -4, path_length) ) h /= np.sqrt(np.sum(h*h)) return h ``` make_path產生一個長度為length，最小延時為delay的指數衰減的波形。這種波形和封閉空間的聲音的傳遞函數有些類似之處，因此在計算機上進行聲音的算法模擬時經常用這種波形作為系統的傳遞函數。: ``` def plot_converge(y, u, label=""): size = len(u) avg_number = 200 e = np.power(y[:size] - u, 2) tmp = e[:int(size/avg_number)*avg_number] tmp.shape = -1, avg_number avg = np.average( tmp, axis=1 ) pl.plot(np.linspace(0, size, len(avg)), 10*np.log10(avg), linewidth=2.0, label=label) def diff_db(h0, h): return 10*np.log10(np.sum((h0-h)*(h0-h)) / np.sum(h0*h0)) ``` plot_converge繪制信號y和信號u之間的誤差，每avg_number個取樣點就上一次誤差的乘方的平均值。我們將用plot_converge函數繪制未知系統的輸出y和自適應濾波器的輸出u之間的誤差。觀察自適應濾波器是如何收斂的，以評價自適應濾波器的收斂特性。diff_db函數同樣是用來評價自適應濾波器的收斂特性，不過他是直接計算未知系統的傳遞函數h0和自適應濾波器的傳遞函數h之間的誤差。下面我們會看到這兩個函數得到的收斂值是相同的。 ### 系統辨識模擬 我們用下面的函數調用nlms算法對系統辨識應用進行模擬： ``` def sim_system_identify(nlms, x, h0, step_size, noise_scale): y = np.convolve(x, h0) d = y + np.random.standard_normal(len(y)) * noise_scale # 添加白色噪聲的外部干擾 h = np.zeros(len(h0), np.float64) # 自適應濾波器的長度和未知系統長度相同，初始值為0 u = nlms( x, d, h, step_size ) return y, u, h ```  系統識別(System Identification)的框圖 此函數的參數分別為： * **nlms** : nlms算法的實現函數 * **x** : 參照信號 * **h0** : 未知系統的傳遞函數，雖然是未知系統，但是計算機模擬時它是已知的 * **step_size** : nlms算法的更新系數 * **noise_scale** : 外部干擾的系數，此系數決定外部干擾的大小，0表示沒有外部干擾 函數的返回值分別為： * **y** : 未知系統的輸出，不包括外部干擾 * **u** : 自適應濾波器的輸出 * **h** : 自適應濾波器的最終的系數 最后我們用下面的函數創建未知系統h0， 參照信號x，然后調用sim_system_identify函數得到結果并且繪圖： ``` def system_identify_test1(): h0 = make_path(32, 256) # 隨機產生一個未知系統的傳遞函數 x = np.random.standard_normal(10000) # 參照信號為白噪聲 y, u, h = sim_system_identify(nlms_numpy.nlms, x, h0, 0.5, 0.1) print diff_db(h0, h) pl.figure( figsize=(8, 6) ) pl.subplot(211) pl.subplots_adjust(hspace=0.4) pl.plot(h0, c="r") pl.plot(h, c="b") pl.title(u"未知系統和收斂后的濾波器的系數比較") pl.subplot(212) plot_converge(y, u) pl.title(u"自適應濾波器收斂特性") pl.xlabel("Iterations (samples)") pl.ylabel("Converge Level (dB)") pl.show() ```  自適應濾波器收斂之后的系數和收斂速度 上部的圖顯示的是未知系統(紅色)和自適應濾波器(藍色)的傳遞函數的系數，我們看到自適應濾波器已經十分接近未知系統了。diff_db(h0, h)的輸出為-25.35dB。下部的圖通過繪制y和u之間的誤差，顯示了自適應濾波器的收斂過程。我們看到經過約3000點的計算之后，收斂過程已經飽和，最終的誤差為-25dB左右，和diff_db計算的結果一致。 從圖中可以看到收斂過程的兩個重要特性：收斂時間和收斂精度。參照信號的特性、外部干擾的大小和更新系數都會影響這兩個特性。下面讓我們看看參照信號為白色噪聲、外部干擾的能量固定時，更新系數對它們影響： ``` def system_identify_test2(): h0 = make_path(32, 256) # 隨機產生一個未知系統的傳遞函數 x = np.random.standard_normal(20000) # 參照信號為白噪聲 pl.figure(figsize=(8,4)) for step_size in np.arange(0.1, 1.0, 0.2): y, u, h = sim_system_identify(nlms_numpy.nlms, x, h0, step_size, 0.1) plot_converge(y, u, label=u"μ=%s" % step_size) pl.title(u"更新系數和收斂特性的關系") pl.xlabel("Iterations (samples)") pl.ylabel("Converge Level (dB)") pl.legend() pl.show() ```  更新系數和收斂速度的關系 下面是更新系數固定，外部干擾能量變化時的收斂特性： ``` def system_identify_test3(): h0 = make_path(32, 256) # 隨機產生一個未知系統的傳遞函數 x = np.random.standard_normal(20000) # 參照信號為白噪聲 pl.figure(figsize=(8,4)) for noise_scale in [0.05, 0.1, 0.2, 0.4, 0.8]: y, u, h = sim_system_identify(nlms_numpy.nlms, x, h0, 0.5, noise_scale) plot_converge(y, u, label=u"noise=%s" % noise_scale) pl.title(u"外部干擾和收斂特性的關系") pl.xlabel("Iterations (samples)") pl.ylabel("Converge Level (dB)") pl.legend() pl.show() ```  外部干擾噪聲和收斂速度的關系 從上面的圖可以看出，當外部干擾的振幅增加一倍、能能量增加6dB時，收斂精度降低6dB。而由于更新系數相同，所以收斂過程中的收斂速度都是一樣的。 ### 信號均衡模擬 對于信號均衡的應用我們用如下的程序進行模擬： ``` def sim_signal_equation(nlms, x, h0, D, step_size, noise_scale): d = x[:-D] x = x[D:] y = np.convolve(x, h0)[:len(x)] h = np.zeros(2*len(h0)+2*D, np.float64) y += np.random.standard_normal(len(y)) * noise_scale u = nlms(y, d, h, step_size) return h ```  信號均衡(Equalization)框圖 sim_signal_equation函數的參數： * **nlms** : nlms算法的實現函數 * **x** : 未知系統的輸入信號 * **h0** : 未知系統的傳遞函數 * **D** : 延遲器的延時參數 * **step_size** : nlms算法的更新系數 * **noise_scale** : 外部干擾的系數，此系數決定外部干擾的大小，0表示沒有外部干擾 在函數中的各個局部變量： * **d** : 輸入信號經過延遲器之后的信號 * **y** : 未知系統的輸出 * **h** : 自適應濾波器的系數，它的長度要足夠長，程序中使用 2倍延時 + 2倍未知系統的傳遞函數的長度 函數的返回值為自適應濾波器收斂后的系數，它能夠均衡h0對輸入信號所造成的影響。我們通過下面的函數產生數據、調用模擬函數以及繪制結果： ``` def signal_equation_test1(): h0 = make_path(5, 64) D = 128 length = 20000 data = np.random.standard_normal(length+D) h = sim_signal_equation(nlms_numpy.nlms, data, h0, D, 0.5, 0.1) pl.figure(figsize=(8,4)) pl.plot(h0, label=u"未知系統") pl.plot(h, label=u"自適應濾波器") pl.plot(np.convolve(h0, h), label=u"二者卷積") pl.title(u"信號均衡演示") pl.legend() w0, H0 = scipy.signal.freqz(h0, worN = 1000) w, H = scipy.signal.freqz(h, worN = 1000) pl.figure(figsize=(8,4)) pl.plot(w0, 20*np.log10(np.abs(H0)), w, 20*np.log10(np.abs(H))) pl.title(u"未知系統和自適應濾波器的振幅特性") pl.xlabel(u"圓頻率") pl.ylabel(u"振幅(dB)") pl.show() ``` 如果延遲器的延時D不夠的話，會由于因果律使得自適應濾波器無法收斂。因此這里我們采用的D的長度為h0的長度的2倍。下圖顯示h0, h和它們的卷積。我們看到h0和h的卷積正好是一個脈沖，其延時為正好等于D(128)。  未知系統和自適應濾波器的級聯(卷積)近似為標準延遲 下圖顯示未知系統的頻率響應(藍色)和自適應濾波器的頻率響應(綠色)，我們看到二者正好相反，也就是說自適應濾波器均衡了未知系統對信號的影響。  未知系統和自適應濾波器的頻率響應正好相反 ### 卷積逆運算 雖然卷積運算最終能歸結為簡單的加法和乘法運算，然而卷積的逆運算就不是很容易計算了。我們知道兩個線性系統h1和h2的級聯h3可以用它們的脈沖響應的卷積計算求得，而所謂卷積的逆運算可以想象為已知h3和h1，求一個h2使它和h1級聯之后正好等于h3。 根據卷積的計算公式可知，如果h1的長度為100，h3的長度為199，那么h2的長度則為100，因為h2的每個系數都是未知的，于是就有100個未知數，而這100個未知數需要滿足199個線性方程：h3中的每個系數都有一個方程與之對應。由于方程數大于未知數的個數，顯然對于任意的h1和h3并不能保證有一個h2使得它和h1的卷積正好等于h3。 既然不能精確求解，那么卷積的逆運算就變成了一個誤差最小化的優化問題。用自適應濾波器計算卷積的逆運算和計算信號均衡類似，將白色噪聲x輸入到h1中得到信號u，將x輸入到h3中得到信號d，然后使用u作為參照信號，d作為目標信號進行NLMS計算，最終收斂后的自適應濾波器的系數就是h2。 下面的程序模擬這一過程： ``` # -*- coding: utf-8 -*- import numpy as np import pylab as pl from nlms_numpy import nlms import scipy.signal as signal def inv_convolve(h1, h3, length): x = np.random.standard_normal(10000) u = signal.lfilter(h1, 1, x) d = signal.lfilter(h3, 1, x) h = np.zeros(length, np.float64) nlms(u, d, h, 0.1) return h h1 = np.fromfile("h1.txt", sep="\n") h1 /= np.max(h1) h3 = np.fromfile("h3.txt", sep="\n") h3 /= np.max(h3) pl.rc('legend', fontsize=10) pl.subplot(411) pl.plot(h3, label="h3") pl.plot(h1, label="h1") pl.legend() pl.gca().set_yticklabels([]) for idx, length in enumerate([128, 256, 512]): pl.subplot(412+idx) h2 = inv_convolve(h1, h3, length) pl.plot(np.convolve(h1, h2)[:len(h3)], label="h1*h2(%s)" % length) pl.legend() pl.gca().set_yticklabels([]) pl.gca().set_xticklabels([]) pl.show() ``` 下面是程序的計算結果：  卷積逆運算演示 程序中的h1和h3從文本文件中讀取而得，它們是ANC(能動噪聲控制)系統中實際測量的脈沖響應。如果能找到一個h2滿足卷積條件的話，就能夠有效的進行噪聲控制。 程序計算出h2的長度分別為128, 256, 512時的結果，可以看出h2越長結果越精確。 ## DLL函數的編寫 ## ctypes的python接口