近期在信號與系統(tǒng)課程中講完了“信號的采樣與恢復(fù)”的內(nèi)容。通常情況下對于信號的采樣都是沿著時間軸對信號的幅值進(jìn)行采樣，獲得信號的離散時間點(diǎn)上的數(shù)據(jù)。

如果將信號 的波形繪制在直接坐標(biāo)系中，那么該曲線就是分布在二維空間上的曲線。曲線上的點(diǎn)可以沿著時間 軸進(jìn)行排列，當(dāng)然也可以按照幅值 的大小進(jìn)行排列。也就是按照取值間隔，將信號通過該間隔的時間進(jìn)行保存，這就是對信號的時間采樣。

沿著時間軸對信號的幅值進(jìn)行采樣，Nyquist-Shannon定理告訴我們?nèi)绾芜M(jìn)行采樣和如何進(jìn)行信號恢復(fù)。但是，如何對信號時間進(jìn)行采樣，如何恢復(fù)是一個經(jīng)典的未解決問題。

Logan定理[Logan,Jr.，1977]對一種特殊信號給出了采樣與恢復(fù)的描述：如果一個信號的頻譜具有倍頻的性質(zhì)，即信號頻譜分布在一個頻率范圍內(nèi)，最高頻率是最低頻率的兩倍。那么這個信號可以通過它的過零點(diǎn)的時間值進(jìn)行恢復(fù)?；謴?fù)的信號與原始信號僅僅相差一個比例因子。

信號采樣與恢復(fù)

1. 幅度采樣

下面是一段普通語音信號的復(fù)制采樣波形。通過簡單的DA轉(zhuǎn)換和低通濾波便可以恢復(fù)出原始的聲音波形。

2. 時間采樣

如果僅僅保留該信號的過零時間點(diǎn)的信息，它的幅值全部去掉，所形成的波形大體上如下圖所示。這很像將原來的語音信號通過一個過零點(diǎn)比較器，輸出的信號反映了信號的極性。

當(dāng)然，這個信號中包含有原來信號的部分信息。但它并不是原來語音信號的完美的回復(fù)。

盡管如此，播放這個信號，還是可以聽到原來語音的信息。雖然有很大的失真和噪聲。這說明原來的信號信息還是部分保留在這些過零點(diǎn)中。

經(jīng)過比較器之后的語音信號為：

如果語音信號滿足Logan定理的要求，那么理論上是可以恢復(fù)出原來的信號的。但如何來恢復(fù)？

如果信號本身是一個周期信號，也就是信號的頻譜是離散的頻譜。相對回復(fù)信號的算法比較簡單。Sam Roweis等人在論文“Signal Reconstruction from Zero-Crossings"中給出了通過求解數(shù)據(jù)矩陣零空間向量的方法，來通過信號的過零時間點(diǎn)來重構(gòu)信號的方法。

重建算法

1. 基本原理

已知到信號 具有倍頻窄帶頻譜，它的頻率范圍分布在 范圍內(nèi)。

已知信號的周期 ， ，以及 個信號過零點(diǎn)：

重構(gòu)信號的計(jì)算步驟如下：

（1）計(jì)算相關(guān)參數(shù)： 以及 ；

（2）構(gòu)造矩陣： ；

（3）尋找數(shù)據(jù)矩陣零空間向量：構(gòu)造數(shù)據(jù)矩陣 ，矩陣大小是 。對該矩陣進(jìn)行奇異值分解，得到 。其中 是 空間上的正交矩陣， 是 空間上的正交矩陣。 是奇異值向量，長度為 。

零空間向量是奇異值向量 中最?。ɡ碚撋蠎?yīng)該為0，但由于計(jì)算誤差的存在，它可能是一個很小的數(shù)）對應(yīng)的 中的向量，由于SVD算法往往把結(jié)果 按照絕對值從大到小排列，所以 的最后一個向量就對應(yīng)著數(shù)據(jù)矩陣的零空間向量。

（4）構(gòu)造信號函數(shù)：將數(shù)據(jù)空間矩陣中零空間向量前 個數(shù)值當(dāng)做 ，后 個數(shù)值當(dāng)做 ，重建信號公式為：

這種回復(fù)信號的過程，實(shí)際上就是根據(jù)信號的過零點(diǎn)來求解上面的函數(shù)中的參數(shù)。具體的理論分析在這里就不再展開了。

2. 測試函數(shù)

（1）實(shí)驗(yàn)信號的數(shù)學(xué)表達(dá)式：

選擇一個頻率分布在10Hz到20Hz之間的一個信號進(jìn)行實(shí)驗(yàn)，隨機(jī)指定對應(yīng)的cos,sin信號的系數(shù)，如下：

這是一個周期為1的倍頻信號。

（2）信號的產(chǎn)生Python程序：

使用下面python程序，可以產(chǎn)生該信號的數(shù)據(jù)。也可以通過該函數(shù)完搜索信號的過零點(diǎn)。

def sfunc1(x): pi2 = 2 * pi retdata = cos(pi2*11*x) + sin(pi2*11*x)/2 + \ cos(pi2*12*x)/33 + sin(pi2*12*x)/4 + \ cos(pi2*13*x) + sin(pi2*13*x)/8 + \ cos(pi2*14*x) + sin(pi2*14*x)/7 + \ cos(pi2*15*x)/22 + sin(pi2*15*x)/3 +\ cos(pi2*16*x) + sin(pi2*16*x)/12 +\ cos(pi2*17*x) + sin(pi2*17*x)/40 +\ cos(pi2*18*x) + sin(pi2*18*x)/2 +\ cos(pi2*19*x)/3 + sin(pi2*19*x)/2
 return retdata

（3）實(shí)驗(yàn)信號的波形：

下面繪制出0~2秒兩個周期內(nèi)的波形。

繪制該信號的過零點(diǎn)飽和信號，它僅僅保留了該信號的過零點(diǎn)的時間和相位信息。計(jì)算公式為：

（4）信號的過零點(diǎn)：

通過數(shù)值計(jì)算，來獲得信號的過零點(diǎn)。下面重新繪制出信號一個周期內(nèi)的波形，沒有添加任何噪聲。

通過對區(qū)間(0,1)采集10^6^個數(shù)值，然后通過尋找過零點(diǎn)，獲得二十八個信號的過零點(diǎn)的值。

搜尋函數(shù)值過零點(diǎn)的python程序如下crosszero(t,val)。其中 是函數(shù)的自變量， 是函數(shù)值的采樣。函數(shù)返回是對應(yīng)函數(shù)過零點(diǎn)時的 的數(shù)值。

def crosszero(t, val): valsign = sign(val) valsignchange = [int(x!=y) for x,y in zip(valsign[0:-1],valsign[1:])] tvalue = [(x,y) for x,y in zip(t[0:-1], valsignchange)] zerot = filter(lambda t: t[1]!= 0, tvalue)
 return [zt[0] for zt in zerot]

通過scipy.optimize.root來尋找信號的根，用于確定信號的過零點(diǎn)。利用上面搜索的結(jié)果作為初始值。

sol = scipy.optimize.root(tssub.sfunc1, czt, method='lm')

• tssub.sfunc1：定義的信號函數(shù)；
• czt：是前面通過數(shù)值過零點(diǎn)搜索獲得的28個根的數(shù)值；

如下是sol['x']中的數(shù)值，包含了最終優(yōu)化后的數(shù)值，對比前面通過搜索獲得數(shù)值，可以看到基本上在10^-6^的內(nèi)存在一定的誤差。

將通過數(shù)值計(jì)算所得到的28個函數(shù)的根繪制在信號波形上，看到他們的分布。

3. 重建結(jié)果

根據(jù)前面所敘述的方法，使用28個過零點(diǎn)信息重構(gòu)出的信號波形如下圖所示。重構(gòu)的信號與原始的信號之間波形基本一致，只是相差了一個比例因子。

前面實(shí)驗(yàn)中的程序和數(shù)據(jù)可以在CSDN博文中看到。

http://zhuoqing.blog.csdn/net