柏路平，馬麗紅，李青龍

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

1 概述

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

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

糾錯(cuò)碼碼字間的最小距離與壓縮采樣矩陣的列相關(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，并且滿(mǎn)足階數(shù)為k ＜p/r +1的RIP 性質(zhì)［8］。代數(shù)曲線矩陣是多項(xiàng)式矩陣的一般形式，選擇合適的代數(shù)曲線能得到比多項(xiàng)式矩陣性能更好的矩陣［9］。二階里德穆勒矩陣［10］將哈達(dá)瑪矩陣作為采樣矩陣的列塊，滿(mǎn)足統(tǒng)計(jì)RIP 性［7］，可利用快速哈達(dá)瑪變換來(lái)加速重構(gòu)。

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

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

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

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

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

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

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

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

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

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

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

3 離散Chirp-Fourier 變換

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6 結(jié)束語(yǔ)

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