突然有個想法，利用機器學習的基本方法——線性迴歸方法，來學習一階RC電路的階躍響應，從而得到RC電路的結構特徵——時間常數τ（即R*C）。回答無疑是肯定的，但問題是怎樣通過最小二乘法、正規方程，以更多的取樣點數來降低訊號採集噪聲對τ估計值的影響。另外，由於最近在搗鼓Jupyter和numpy這些東西，正好嘗試不用matlab而用Jupyter試試看。結果是意外的好用，尤其是在Jupyter指令碼中插入LaTeX格式的公式的功能，真是太方便了！嘗試了直接把紙上手寫的公式轉換到Jupyter指令碼中的常見工具軟體。

以下原創內容歡迎網友轉載，但請註明出處： https://www.cnblogs.com/helesheng

一、理論推導

1．線性迴歸分析及正規方程

傳統意義說，線性迴歸問題是用最小二乘法（即正規方程），解決線性方程組的均方誤差最小化問題。已知輸出輸入X是由多個變數構成的行向量，W是係數向量（列向量），b為偏置

 在機器學習中，把每次的輸入x作為一行組成更大的矩陣，即每一行代表一個樣本，該矩陣稱為設計矩陣X^(train)。若樣本數為k，則X^(train)有k行，每個樣本都會得到一個輸出y，將y集合成一個列向量Y^(train)，k個相同的b也組成列向量b。為簡化表達，將b簡化為附加在係數列向量W最後的常數b，X^(train)則在每行的最後增加一個1，W（列向量）的最後增加一個待估計的b。為了使估計的結果：

  和Y^(train)之差的平方和最小，有正規方程可以求解W:

 2.一階RC電路的階躍響應

一階RC電路的電路模型如下圖所示。

圖1 一階RC電路

幅度為Vcc的階躍訊號從Vin處輸入，在Vout處測量輸出。解微分方程可得自變數為時間t的響應函式。 

 其中時間常數τ = R*C。我希望通過測量階躍訊號輸入條件下，實際RC電路的響應曲線V(t)，並通過V(t)來估計時間常數τ。如果做純理論推導，只要知道任意時刻t0的輸出電壓V(t0)就可以解方程（2）得到τ。但在實際電路中電壓V(t0)的測量往往受到諸多幹擾的影響，並不準確。是否可以測量多個t值時刻對應的V(t)，並利用正規方程（1）得到一個統計意義上最優的估計呢？是接下來要解決的問題。

3.非線性函式的最小二乘估計

仔細觀察適用正規方程的目標函式（0）式的特點，可以發想讓非線性的要讓（2）式能夠使用正規方程，必須讓：

1）     含有待估計的變數τ的函式充當（0）式中的係數W，而設計矩陣X^(train)則可以由含有時間t或測量電壓V(t)的函式充當。

2）     W和X^(train)之間的關係必須是簡單的相乘。

 顯然，只有用時間t的序列充當設計矩陣X^(train)，才有可能使W和X^(train)之間的關係必須是相乘。至於其他的非線性部分則可以通過等效變換，變換到等式的另一側來充當Y^(train)。綜上，可以將（2）式變換為（3）式。

（3）式的整個左邊可以計算得到Y^(train)；（3）式右邊的時間t的序列在構成設計矩陣X^(train)，1/τ則構成相當於（0）式中的係數矩陣W。這樣就可以通過正規方程（2）式來求解統計意義上最優的估計了。當然，從（3）式也可以看出，經過線性校正的模型中係數向量W只有一個元素——是個標量，偏置b也是恆等於0的。

 二、模擬模型

想利用最近正在嘗試使用的Jupyter和numpy替代熟悉的Matlab，驗證上述非線性函式最小二乘估計的想法。下面先建立一個模型：

1）     輸入為幅度Vcc為3.3V的階躍訊號；

2）     時間常數τ為0.1秒；

3）     簡單模擬取樣間隔的隨機性：是間隔等於delta1=0.0015秒和delta1=0.0011秒的兩個序列的疊加。

4）     取樣總長度為n=5倍τ；

5）     訊號上疊加了幅度約為20mV的白噪聲——至於為什麼是20mV，將在後續部分詳細介紹。

利用python和numpy進行數值模擬的程式碼如下：

 1 import numpy as np
 2 import matplotlib.pyplot as plt
 3 tao=0.1#時間常數
 4 vcc=3.3#電源電壓
 5 n=5#時長取時間常數tao的n倍
 6 delta1=0.0015#第一種取樣間隔
 7 delta2=0.0011#第一種取樣間隔
 8 t1=delta1*np.arange(np.ceil(n*tao/delta1))
 9 t2=delta2*np.arange(np.ceil(n*tao/delta2))
10 t=np.append(t1,t2)#為演示最小二乘擬合的功能，故意結合兩種取樣率下的時間點
11 t.sort()#對t進行排序
12 plt.plot(t)
13 s=vcc*(1-np.exp(-t/tao))#理論的波形曲線
14 plt.plot(t,s)#注意這裡的plot函式使用了x軸和y軸兩個軸，因為s中的資料不是均勻的
15 N_amp=np.exp(-n)*vcc
16 N_amp
17 noise = np.random.uniform(-N_amp, N_amp, (len(t)))#噪聲：正負np.exp(-5)*3.3之間均勻分佈
18 s_nr=s+noise#加入噪聲後的訊號
19 plt.plot(t,s_nr)
20 yt=np.log(vcc/(vcc-s_nr))
21 plt.plot(t,yt)
22 yt=np.mat(yt)
23 yt=yt.T
24 x=np.mat(t)#將X轉換為矩陣資料型別
25 x=x.T#正規方程中的x應該是個列向量
26 w=(np.linalg.inv(x.T*x))*x.T*yt#求解正規方程
27 E_tao = np.round(1/float(w),4)#對時間常數的tao的估計，保留到4位小數
28 E_tao

說明：

1）     其間序列包含了兩個等差數列t1和t2的融合，它們的間隔互質，沒有重複，目的是模擬取樣時間的隨機性。最後用sort()方法又對時間序列進行排序的目的是為了後續容易觀察波形更直觀。如果僅僅為了使用正規方程，是不需要重新排序的。

2）     校正的非線性方程（3）和原始方程（2）相比有一個重大的缺陷就是：左側求對數的括號內的數值不能為負——如果是純理論推導，這是不可能發生的。但假如隨機噪聲後的V(t)是有可能大於階躍幅度Vcc的，此時括號內將變為一個負數，使得計算變得沒有意義。我的解決之道是將假如的隨機噪聲幅度限制在模擬截止時刻V(t)和Vcc之差的範圍內，及程式碼中的N_amp。由於模擬的結束時刻為n（=5）個τ，所以N_amp的值等於np.exp(-n)*vcc。
這樣做沒有任何理論依據，僅僅是受限於（3）和（2）式之間的非完全等價變換——屬於線性化校正需要付出的代價。

3）     由於待估計的引數只有一個（1/τ）所以正規方程的解也是隻有一個元素的矩陣。將其轉換為標量後取倒得到最優估計。

三、結果評估

加入噪聲後的訊號如下圖所示，與通常情況的實測波形吻合度很高。

圖2 模擬產生的帶有噪聲的階躍響應 

 對這些加入噪聲的訊號進行線性校正後得到進行最小二乘估計的訊號yt為下圖所示。其基本趨勢是一條斜率為（1/τ）的直線，和我預計的結果一樣。

圖3 對圖2進行線性校正後的待估計訊號

但可以明顯的看到，由於（3）式左側在V(t)的大小接近Vcc時對加入的白噪聲進行了放大。因此當t遞增時，由白噪聲造成的訊號的不確定性大大增加了。也就是在套用正規方程時，t值較大時的噪聲對結果的影響大於t值較小時的噪聲對結果的影響。這是使用非線性校正函式需要付出的重要代價。

另外，多次執行以上程式碼的得到 都是一個略小於實際τ（=0.1）的數值——約為0.099左右，也就是說， 不是無偏估計。這應該是由於線性校正函式（（3）式左側），對於噪聲noise大於0和小於0的部分進行了非對稱的變換造成的。這雖然造成的偏差不大，但也是使用非線性校正函式需要付出的代價。 

四、Jupyter notebook

上述練習的一個重要目的是練習使用Jupyter notebook，並在其中內嵌具有很好互動性的公式等資訊。以下是部分程式執行效果的截圖，雖然我對markdown語法還不熟悉，格式和排版還不夠漂亮，但還是能夠明顯的提高程式碼的可讀性。

圖4 Jupyter notebook執行效果截圖

 其中需要重點記錄下的是公式程式碼的嵌入過程：

1）我首先在紙上手寫了一些公式，用手機對其拍照，如： 

圖5 手寫的公式

  2）用mathpix tools對這些照片截圖，並掃描（mathpix tools有windows版和iOS版，均可免費試用）。Mathpix可以直接得到LaTeX格式的公式表示式。下圖是iOS版本的mathpix介面截圖。

圖6 iOS版本的mathpix截圖

3）在Jupyter notebook上部的下拉選單中選擇單元格的格式為Markdown；將LaTeX格式的公式表示式貼上到該單元格內，並在LaTeX公式表示式的前後加上“$$”符號，執行該單元格即可得到圖4所示效果的公式了。

圖7 在Jupyter notebook中輸入LaTeX公式

 五、進一步的實際測試

工作做到這裡離其實並沒有完，還應該做一個簡單的實際電路實測一下。我會在後續的假期中抽半天時間，在STM32開發板上完成這個工作：用GPIO產生一個節約訊號後，連續採集5個τ時間長度的訊號，並代入正規方程求解，歡迎大家繼續關注更新。

…