利用機器學習的基本方法——線性回歸方法，來學習一階RC電路的階躍響應，進而得到RC電路的結構特征——時間常數τ（即R*C）。回答無疑是肯定的，但問題是怎樣通過最小二乘法、正規方程，以更多的采樣點數來降低信号采集噪聲對τ估計值的影響。另外，由于最近在搗鼓Jupyter和numpy這些東西，正好嘗試不用matlab而用Jupyter試試看。結果是意外的好用，尤其是在Jupyter腳本中插入LaTeX格式的公式的功能，真是太友善了！嘗試了直接把紙上手寫的公式轉換到Jupyter腳本中的常見工具軟體。
突然有個想法，利用機器學習的基本方法——線性回歸，來學習一階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個τ時間長度的信号，并代入正規方程求解，歡迎大家繼續關注更新。
……