近期在信號與系統課程中講完了“信號的採樣與恢復”的內容。通常情況下對於信號的採樣都是沿著時間軸對信號的幅值進行採樣，獲得信號的離散時間點上的數據。
如果將信號的波形繪製在直接座標系中，那麼該曲線就是分佈在二維空間上的曲線。曲線上的點可以沿著時間軸進行排列，當然也可以按照幅值的大小進行排列。也就是按照取值間隔，將信號通過該間隔的時間進行保存，這就是對信號的時間採樣。
聲音的時間波形
沿著時間軸對信號的幅值進行採樣，Nyquist-Shannon 定理告訴我們如何進行採樣和如何進行信號恢復。但是，如何對信號時間進行採樣，如何恢復是一個經典的未解決問題。
Logan 定理[Logan,Jr.，1977]對一種特殊信號給出了採樣與恢復的描述：如果一個信號的頻譜具有倍頻的性質，即信號頻譜分佈在一個頻率範圍內，最高頻率是最低頻率的兩倍。那麼這個信號可以通過它的過零點的時間值進行恢復。恢復的信號與原始信號僅僅相差一個比例因子。
信號採樣與恢復
1. 幅度採樣
下面是一段普通語音信號的複製採樣波形。通過簡單的 DA 轉換和低通濾波便可以恢復出原始的聲音波形。
幅度採集聲音波形
這段語音信號的內容為：
2. 時間採樣
如果僅僅保留該信號的過零時間點的信息，它的幅值全部去掉，所形成的波形大體上如下圖所示。這很像將原來的語音信號通過一個過零點比較器，輸出的信號反映了信號的極性。
<code>聲音信號過零點採樣/<code>
當然，這個信號中包含有原來信號的部分信息。但它並不是原來語音信號的完美的回覆。
儘管如此，播放這個信號，還是可以聽到原來語音的信息。雖然有很大的失真和噪聲。這說明原來的信號信息還是部分保留在這些過零點中。
經過比較器之後的語音信號為：
如果語音信號滿足 Logan 定理的要求，那麼理論上是可以恢復出原來的信號的。但如何來恢復？
如果信號本身是一個週期信號，也就是信號的頻譜是離散的頻譜。相對回覆信號的算法比較簡單。Sam Roweis 等人在論文“Signal Reconstruction from Zero-Crossings"中給出了通過求解數據矩陣零空間向量的方法，來通過信號的過零時間點來重構信號的方法。
重建算法
1. 基本原理
已知到信號具有倍頻窄帶頻譜，它的頻率範圍分佈在範圍內。
倍頻信號頻譜示意圖
已知信號的週期，，以及個信號過零點：
重構信號的計算步驟如下：
（1）計算相關參數： 以及 ；
（2）構造矩陣： ；
（3）尋找數據矩陣零空間向量：構造數據矩陣，矩陣大小是。對該矩陣進行奇異值分解，得到。其中是空間上的正交矩陣，是空間上的正交矩陣。是奇異值向量，長度為。
零空間向量是奇異值向量中最小（理論上應該為 0，但由於計算誤差的存在，它可能是一個很小的數）對應的中的向量，由於 SVD 算法往往把結果按照絕對值從大到小排列，所以的最後一個向量就對應著數據矩陣的零空間向量。
（4）構造信號函數：將數據空間矩陣中零空間向量前個數值當做，後個數值當做，重建信號公式為：
這種回覆信號的過程，實際上就是根據信號的過零點來求解上面的函數中的參數。具體的理論分析在這裡就不再展開了。
2. 測試函數
（1）實驗信號的數學表達式：
選擇一個頻率分佈在 10Hz 到 20Hz 之間的一個信號進行實驗，隨機指定對應的 cos,sin 信號的係數，如下：
信號的數學表達式
這是一個週期為 1 的倍頻信號。
（2）信號的產生 Python 程序：
使用下面 python 程序，可以產生該信號的數據。也可以通過該函數完搜索信號的過零點。
<code>def sfunc1(x):/<code><code> pi2 = 2 * pi/<code><code> retdata = cos(pi2*11*x) + sin(pi2*11*x)/2 + \/<code><code> cos(pi2*12*x)/33 + sin(pi2*12*x)/4 + \/<code><code> cos(pi2*13*x) + sin(pi2*13*x)/8 + \/<code><code> cos(pi2*14*x) + sin(pi2*14*x)/7 + \/<code><code> cos(pi2*15*x)/22 + sin(pi2*15*x)/3 +\/<code><code> cos(pi2*16*x) + sin(pi2*16*x)/12 +\/<code><code> cos(pi2*17*x) + sin(pi2*17*x)/40 +\/<code><code> cos(pi2*18*x) + sin(pi2*18*x)/2 +\/<code><code> cos(pi2*19*x)/3 + sin(pi2*19*x)/2/<code><code> return retdata/<code>
（3）實驗信號的波形：
下面繪製出 0~2 秒兩個週期內的波形。
測試函數 sfunc1 信號波形
繪製該信號的過零點飽和信號，它僅僅保留了該信號的過零點的時間和相位信息。計算公式為：
該信號的幅值飽和信號
（4）信號的過零點：
通過數值計算，來獲得信號的過零點。下面重新繪製出信號一個週期內的波形，沒有添加任何噪聲。
一個週期（0~1）之間的信號波形
通過對區間(0,1)採集 10^6^個數值，然後通過尋找過零點，獲得二十八個信號的過零點的值。
搜尋函數值過零點的 python 程序如下 crosszero(t,val)。其中是函數的自變量，是函數值的採樣。函數返回是對應函數過零點時的的數值。
<code>def crosszero(t, val):/<code><code> valsign = sign(val)/<code><code> valsignchange = [int(x!=y) for x,y in zip(valsign[0:-1],valsign[1:])]/<code><code> tvalue = [(x,y) for x,y in zip(t[0:-1], valsignchange)]/<code><code> zerot = filter(lambda t: t[1]!= 0, tvalue)/<code><code> return [zt[0] for zt in zerot]/<code>
通過 scipy.optimize.root 來尋找信號的根，用於確定信號的過零點。利用上面搜索的結果作為初始值。
<code>sol = scipy.optimize.root(tssub.sfunc1, czt, method='lm')/<code>
- tssub.sfunc1：定義的信號函數；czt：是前面通過數值過零點搜索獲得的 28 個根的數值；
如下是 sol['x']中的數值，包含了最終優化後的數值，對比前面通過搜索獲得數值，可以看到基本上在 10^-6^的內存在一定的誤差。
將通過數值計算所得到的 28 個函數的根繪製在信號波形上，看到他們的分佈。
尋找到的信號過零點
3. 重建結果
根據前面所敘述的方法，使用 28 個過零點信息重構出的信號波形如下圖所示。重構的信號與原始的信號之間波形基本一致，只是相差了一個比例因子。
重建的波形結果
前面實驗中的程序和數據可以在 CSDN 博文中看到。
http://zhuoqing.blog.csdn/net
