柏路平，馬麗紅，李青龍

(華南理工大學(xué)電子與信息學(xué)院，廣州 510640)

1 概述

壓縮采樣的目的主要是用盡可能少的采樣數(shù)據(jù)獲取稀疏信號［1-2］，它將一個k-稀疏(非零元個數(shù)為k)信號向量X 在采樣矩陣Φ 上投影，得到一個只有少量數(shù)據(jù)的采樣信號向量Y?？蓧嚎s非稀疏信號在某些正交變換下只保留少量的大系數(shù)，也可視為稀疏信號，其余去掉的接近于零的小系數(shù)，可看作信號上疊加的噪聲。

從X 得到Y(jié) 是一個簡單的線性投影過程，然而從Y 準(zhǔn)確重構(gòu)X 卻是一個通過最佳原子選擇估計(jì)信號的非線性過程。首先，最小l1范數(shù)優(yōu)化［3］和l0范數(shù)貪婪迭代算法［4］是常用的重構(gòu)方法;其次，只有任意列之間的相關(guān)系數(shù)很小的矩陣才能成為采樣矩陣，并且需要滿足k 階受限等距性質(zhì)(Restricted Isometry Property，RIP)［5］，才能以較大的概率保證準(zhǔn)確重構(gòu)X。RIP 性質(zhì)要求Φ 的任意M 列組成的子矩陣具有近似正交矩陣的性質(zhì)。在RIP 性質(zhì)上隨機(jī)矩陣是Φ 的最佳選擇，高斯矩陣和伯努力矩陣列之間的相關(guān)性很小，能以最小采樣長度滿足RIP 性質(zhì)［6］，但其在硬件電路中難以實(shí)現(xiàn)，需要很大的存儲空間存放采樣矩陣的隨機(jī)元素，以及需要大量重復(fù)實(shí)驗(yàn)驗(yàn)證其統(tǒng)計(jì)性能［7］。

糾錯碼碼字間的最小距離與壓縮采樣矩陣的列相關(guān)系數(shù)相似，因此，線性碼成為確定性采樣矩陣的主要構(gòu)造方法。多項(xiàng)式矩陣是素?cái)?shù)域GF(p)的r 階多項(xiàng)式構(gòu)造p2×pr+1的二進(jìn)制稀疏矩陣，其任意兩列相關(guān)系數(shù)最大為r/p，并且滿足階數(shù)為k ＜p/r +1的RIP 性質(zhì)［8］。代數(shù)曲線矩陣是多項(xiàng)式矩陣的一般形式，選擇合適的代數(shù)曲線能得到比多項(xiàng)式矩陣性能更好的矩陣［9］。二階里德穆勒矩陣［10］將哈達(dá)瑪矩陣作為采樣矩陣的列塊，滿足統(tǒng)計(jì)RIP 性［7］，可利用快速哈達(dá)瑪變換來加速重構(gòu)。

由離散Chirp 序列構(gòu)成的Chirp 矩陣采樣，重構(gòu)時可利用離散傅里葉變換快速算法(FFT)［11］。但由于chirp 信號交調(diào)干擾的存在，用離散傅立葉變換(DFT)相關(guān)檢測方法，即用采樣信號自相關(guān)的DFT最大的幅值及位置來估計(jì)Chirp 參數(shù)，Chirp 率、基頻和幅值的估計(jì)誤差會隨信號稀疏度k 的增大而急劇增加。

本文提出基于DCFT 的Chirp 矩陣采樣重構(gòu)算法，在保留Chirp 矩陣的確定性，以及快速傅立葉變換加速重構(gòu)的優(yōu)點(diǎn)同時，能夠增大可采樣信號的稀疏度k 以及提高信號的重構(gòu)精度。DCFT 將Chirp信號的頻譜從一維擴(kuò)展至二維，使Chirp 率相同而基頻不同的分量不再相互干擾，比DFT 具有更高的頻譜分辨率［12］。用DCFT 替代DFT 相關(guān)法來檢測采樣信號的最佳原子，可避免DFT 相關(guān)檢測算法因交調(diào)干擾造成的原子誤檢測，最后用最小二乘法估計(jì)非零元的幅值。

2 DFT 相關(guān)檢測重構(gòu)算法的局限性

基于線性調(diào)頻原理的Chirp 采樣矩陣［11］，是由Chirp 率和基頻順序增加的離散Chirp 序列排列而成。設(shè)c，f 分別表示Chirp 率和基頻，長為M 的離散Chirp 序列向量Vc，f可表示為:

其中，(c，f)在整數(shù)域?M中共有M2種組合，將這M2個Chirp 序列按式(2)的順序構(gòu)成M × M2矩陣Φ:

如果X 是k 稀疏的，則經(jīng)矩陣Φ 采樣得到的信號Y=ΦX 可表示為k 個chirp 序列的線性組合公式如下:

其中，Xi表示X 的第i 個非零元，則采樣信號Y 的互相關(guān)向量F 具有下面的形式:

其中，時移T∈?M且T≠0，交叉項(xiàng)具有形式:

共有k(k-1)項(xiàng)，將引入交調(diào)干擾，其頻譜不規(guī)則地分布在所有DFT 系數(shù)上。

DFT 相關(guān)檢測重構(gòu)算法的原理是:互相關(guān)向量F 的頻譜在2ciTmod(M)的離散頻率處為尖脈沖，如果M 是素?cái)?shù)，那么從chirp 率到F 的DFT 系數(shù)是一個一一映射。所以，對F 做DFT 變換，在2ciTmod(M)的位置有一個峰值，可得到Chirp 率ci。然后將Y 乘以e-j2πcil2/M2，去除Chirp 率ci，Y 就變成具有單一諧波fi分量的信號，其DFT 峰值位置和大小分別對應(yīng)一個基頻fi和幅值Xci，fi，由此得到一個非零元的位置索引估計(jì)Mci+fi和幅值估計(jì)Xci，fi。然后將Y 用殘差Y-Xci，fiej2πcil2/Mej2πfil/M代替，經(jīng)過k 次迭代，就可按照幅值從大到小逐步得到Xi的位置與幅值。根據(jù)式(4)的時移T 的不同取值，文獻(xiàn)［11］給出單時移和多時移的DFT 相關(guān)檢測2 種重構(gòu)算法，多時移DFT 相關(guān)重構(gòu)算法具有更小的重構(gòu)誤差以及更強(qiáng)的抗噪性能。

用Chirp 矩陣采樣時，DFT 相關(guān)檢測算法重構(gòu)的信號誤差很大，并且對稀疏度k 較大的信號會因誤差太大而無法重構(gòu)。(1)DFT 算法的采樣數(shù)為信號長度的根方可采樣的信號稀疏度k 要比采樣數(shù)小;(2)chirp 信號交調(diào)干擾的存在，DFT 相關(guān)檢測算法沒有考慮式(4)中的交叉項(xiàng)的影響。隨著信號稀疏度k 的增加，組成采樣信號的Chirp 分量個數(shù)也相應(yīng)增加，而式(5)那樣的交叉項(xiàng)個數(shù)卻成平方增加，對重構(gòu)結(jié)果的影響急劇增加。再用DFT 相關(guān)算法檢測最佳原子以及簡單地用頻域fi幅值的平方根來估計(jì)非零元的幅值，重構(gòu)誤差也將急劇增加。

3 離散Chirp-Fourier 變換

DFT 將時域的離散諧波分量無交叉干擾地映射到頻域，每個離散頻率都與一個諧波一一對應(yīng)，并且存在FFT 快速算法，所以DFT 成為頻譜估計(jì)的主要手段。然而DFT 對于頻率線性時變的Chirp 信號的參數(shù)估計(jì)卻不太理想，因?yàn)椴煌珻hirp 率和基頻的分量在頻域都有重疊。DCFT 將DFT 從一維擴(kuò)展為二維，是DFT 的一般形式:

其中，M 為DCFT 變換的長度，WM=e-j2π/M;c 表示Chirp 率;f 表示基頻。容易知道，當(dāng)c 固定時，{X(c，f)}0≤f≤M -1為X(n)Wcn2M 的DFT 變換，因此FFT 也可用于X(c，f)的計(jì)算。DCFT 的一個重要性質(zhì)為:

圖1 (1×n2+3×n)長度為5 的DCFT 幅值分布

4 Chirp 矩陣采樣重構(gòu)算法

由式(4)和式(5)可以看出，采樣信號Y 的所有Chirp 分量在其互相關(guān)向量F 的DFT 域都有交調(diào)干擾，而Chirp 率相同而基頻不同的分量在DCFT 域中沒有疊加，所以用DCFT 替代DFT 相關(guān)算法可以更準(zhǔn)確地檢測采樣信號的最佳原子，進(jìn)而減小重構(gòu)誤差。另外，在壓縮采樣中，稀疏度k 是信號的重要信息，為保證準(zhǔn)確重構(gòu)，采樣數(shù)應(yīng)隨k 的增大而增大。在DFT 相關(guān)檢測的重構(gòu)算法中，采樣數(shù)由信號長度N 決定，無論稀疏度為何值，采樣數(shù)均為可采樣的信號的稀疏k 非常有限。在這種情況下，增加采樣數(shù)，以降低壓縮效率換取重構(gòu)性能是準(zhǔn)確重構(gòu)稀疏度大的信號的主要方法。

當(dāng)M 為素?cái)?shù)時，在單分量Chirp 信號的M 長DCFT 中，坐標(biāo)(c0，f0)處的幅值為非(c0，f0)處最大為1。假設(shè)信號由等幅值的k 個Chirp 分量組成，當(dāng)滿足時，|X(c，f)|中最大k 個幅值就能以很大的概率分別對應(yīng)所有的k 個Chirp 分量。用DCFT 來估計(jì)采樣信號Y 的Chirp 參數(shù)時，變換的長度與采樣數(shù)相等，所以為利用Y 的DCFT 變換中k 個最大的幅值所對應(yīng)的原子索引來擊中非零元的位置，采樣數(shù)M 應(yīng)增大為k2，但只有當(dāng)作DCFT 變換的長度為素?cái)?shù)時，才具有式(8)的性質(zhì)，所以最終的采樣數(shù)為大于k2的最小素?cái)?shù)。

根據(jù)前面的分析，本文基于DCFT 的Chirp 矩陣采樣重構(gòu)算法的采樣矩陣Φ 的元素為:

DCFT 重構(gòu)算法的執(zhí)行步驟為:

(1)計(jì)算采樣信號Y 長為M 的DCFT 變換X(c，f);

(2)搜索|X(c，f)|中最大的k 個值的坐標(biāo)(ci，fi)，記錄位置索引集合Λ=∪k-1i=0{M·ci+fi}，從采樣矩陣Φ 中提取Λ 對應(yīng)的原子組成最佳匹配矩陣ΦΛ;

(3)用最小二乘法估計(jì)重構(gòu)信號的非零元向量XΛ的幅值:

(4)在長度為N 的全零向量中位置索引為Λ 中的元素依次取XΛ中的元素就可得到信號X 的估計(jì)值Xt。

與DFT 相關(guān)檢測重構(gòu)算法相比，本文基于DCFT 的重構(gòu)算法主要有3 個方面的改進(jìn):(1)用頻譜分辨率更高的DCFT 代替DFT 相關(guān)算法檢測采樣信號的最佳原子和非零元的位置;(2)適當(dāng)增加采樣數(shù)至信號稀疏度的平方以換取重構(gòu)性能;(3)提取最佳原子匹配矩陣ΦΛ，用最小二乘法求解非零元的幅度值，使重構(gòu)的Xt與X 之間的誤差達(dá)到最小。

5 實(shí)驗(yàn)性能與分析

為了驗(yàn)證DCFT 重構(gòu)算法的性能，本文用Matlab 程序仿真一維稀疏信號的采樣與重構(gòu)過程，并用重構(gòu)信號Xt與X 的相對誤差:

來評估重構(gòu)誤差。為了與DFT 重構(gòu)算法比較對相同長度和稀疏度的信號的重構(gòu)性能，選擇長度N=412=1 681 的測試信號，對于每個稀疏度k ∈{1，2，…，40}，非零元位置隨機(jī)產(chǎn)生而幅值服從(0，1)上均勻分布的信號，用本文提出的采樣和重構(gòu)方法求得重構(gòu)誤差SE，并重復(fù)1 000 次實(shí)驗(yàn)，得到均方誤差MSE。DCFT 重構(gòu)算法的均方誤差—稀疏度性能曲線如圖2 所示。

圖2 DCFT 重構(gòu)算法的誤差曲線

為便于比較，2 種DFT 相關(guān)檢測重構(gòu)算法以及它們用最小二乘法優(yōu)化后的誤差性能曲線也在圖中畫出。優(yōu)化過程是先用DFT 相關(guān)算法檢測最佳原子，再用最小二乘法求解非零元的幅值，這樣重構(gòu)的信號具有更小的重構(gòu)誤差。

在圖2 中，SS 和MS 分別表示單時移與多時移DFT 相關(guān)檢測算法，SSLS 和MSLS 分別表示用最小二乘法優(yōu)化后的單時移與多時移DFT 相關(guān)檢測算法的誤差曲線?？梢钥闯?，在稀疏度k 分別大于6和12 時，最小二乘優(yōu)化后的單時移和多時移的DFT相關(guān)檢測算法的相對均方誤差已大于0.1，可視為重構(gòu)完全失效。而DCFT 重構(gòu)算法只在k=6 時有誤差約為0.09 峰值，之后隨k 的增加，誤差逐漸減小。如果將誤差小于0.1 的重構(gòu)視為準(zhǔn)確重構(gòu)，則DCFT重構(gòu)算法能準(zhǔn)確重構(gòu)的信號稀疏度k 擴(kuò)大為DFT 算法的4 倍左右。

由于Chirp 采樣矩陣可由信號長度N 和稀疏度k 唯一確定，重構(gòu)時不需存儲空間大小為M ×N 的采樣矩陣，只需存儲長為N 的DCFT 變換矩陣及重構(gòu)信號Xt，因此重構(gòu)空間復(fù)雜度為O(N)。FFT 變換的時間復(fù)雜度為O(MlogM)，而計(jì)算DCFT 變換只需要次FFT 變換，所以復(fù)雜度為O(Nlogk)，從N 個數(shù)里面搜索最大的k 個數(shù)的排序時間復(fù)雜度為O(kN)，解最小二乘的復(fù)雜度為O(k2M)，所以整個DCFT 重構(gòu)算法的時間復(fù)雜度為O(kN)。當(dāng)信號長為O時，其復(fù)雜度處于單時移DFT 相關(guān)檢測的O(kMlogM)和多時移DFT 相關(guān)檢測的O(kM2logM)之間，但優(yōu)于最小l1范數(shù)中內(nèi)點(diǎn)法O(M2N3/2)以及匹配追蹤的O(kMN)時間復(fù)雜度。

6 結(jié)束語

本文提出基于DCFT 的Chirp 矩陣采樣重構(gòu)算法，采樣時根據(jù)信號X 的稀疏度k，增加采樣數(shù)以換取重構(gòu)信號的質(zhì)量;重構(gòu)時，搜索采樣信號DCFT 變換最大的k 個幅值，其坐標(biāo)所對應(yīng)的在采樣矩陣中的原子即為信號X 的最佳匹配原子，最后用最小二乘法求解非零元素的幅值。對一維信號的采樣和重構(gòu)實(shí)驗(yàn)結(jié)果表明，與DFT 相關(guān)檢測算法相比，DCFT重構(gòu)算法具有更小重構(gòu)誤差。本文提出方法的采樣數(shù)隨信號稀疏度呈平方指數(shù)增加，雖然能準(zhǔn)確重構(gòu)稀疏度k 更大的信號，但與最佳的線性采樣數(shù)相差較遠(yuǎn)，因此，保證重構(gòu)質(zhì)量的同時減少采樣數(shù)成為下一步的工作目標(biāo)。

［1］Duarte M F，Eldar Y C.Structured Compressed Sensing From Theory to Applications［J］.IEEE Transactions on Signal Processing，2011，59(9):4053-4085.

［2］馬堅(jiān)偉，徐 杰，鮑躍全，等.壓縮感知及其應(yīng)用:從稀疏約束到低秩約束優(yōu)化［J］.信號處理，2012，28(5):609-623.

［3］Becker S，Bobin J，Candes E.NESTA:A Fast and Accurate First-order Method for Sparse Recovery［J］.SIAM Journal on Imaging Sciences，2011，4(1):1-39.

［4］Dai W，Milenkovic O.Subspace Pursuit for Compressive Sensing Signal Reconstruction［J］.IEEE Transactions on Information Theory，2009，55(5):2230-2249.

［5］Davenport M A，Wakin M B.Analysis of Orthogonal Matching Pursuit Using the Restricted Isometry Property［J］.IEEE Transactions on Information Theory，2010，56(9):4395-4401.

［6］Candes E J，Tao T.Near-optimal Signal Recovery Form Random Projections:Universal Encoding Strategies［J］.IEEE Transactions on Information Theory，2006，52(12):5406-5425.

［7］Calderbank R，Howard S，Jafarpour S.Construction of a Large Class of Deterministic Sensing Matrices That Satisfy a Statistical Isometry Property［J］.IEEE Journal of Selected Topics in Signal Processing，2010，4(2):358-374.

［8］Devore R A.Deterministic Constructions of Compressed Sensing Matrices［J］.Journal of Complexity，2007，23(4-6):918-925.

［9］Li Shuxing，Gao Fei，Ge Gennian.Deterministic Construction of Compressed Sensing Matrices via Algebraic Curves［J］.IEEE Transactions on Information Theory，2012，58(8):5035-5041.

［10］Howard S D，Calderbank A R，Searle S J.A Fast Reconstruction Algorithm for Deterministic Compressive Sensing Using Second Order Reed-Muller Codes［C］//Proceedings of the 42nd Annual Conference on Information Sciences and Systems.Princeton，USA:IEEE Press，2008:231-239.

［11］Applebaum L，Howard S，Searle S.ChirpSensing Codes:Deterministic Compressed Sensing Measurements for Fast Recovery［J］.Applied and Computational Harmonic Analysis，2009，26(1):283-290.

［12］Xia Xianggen.Discrete Chirp-fourier Transform and Its Application to Chirp Rate Estimation ［J］.IEEE Transactions on Signal Processing，2000，48 (11):3122-3133.