突然有個想法，利用機器學習的基本方法——線性回歸，來學習一階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個τ時間長度的信號，並代入正規方程求解，歡迎大家繼續關注更新。
……