基本概念
關(guān)于壓縮感知的一些基礎(chǔ)概念倔毙，可以看 Baraniuk的這篇Lecture Note（想不到都有4000+的引用量）果元，還有Baraniuk在Microsoft Research的一個(gè)報(bào)告坐漏，以及一個(gè)基本和報(bào)告內(nèi)容一致的PDF文件
Shannon/Nyquist采樣定律
Shannon/Nyquist采樣定律 告訴我們，如果要在采樣以后能夠完整保留保證原始信號(hào)中信息，采樣頻率必須大于信號(hào)中最高頻率的2倍。一直以來(lái)粒氧，聲音、圖像等傳感器采樣的方式都服從了Shannon/Nyquist采樣定律节腐，在保證了信號(hào)完整性的同時(shí)帶來(lái)的是大量的數(shù)據(jù)外盯。比如圖像傳感器采集的原始圖像數(shù)據(jù)，每張圖都有好幾十兆铜跑。為了減小數(shù)據(jù)量门怪，才出現(xiàn)了諸如JPEG，JPEG2000等壓縮算法锅纺。
圖像來(lái)源
這樣的過(guò)程可以用上圖的上半部分來(lái)表示，x為原始信號(hào)肋殴，經(jīng)過(guò)Uniform sample 采樣得到了一個(gè)N維度的向量囤锉，然后在經(jīng)過(guò)DCT或者小波變換，得到了M個(gè)主要的成分护锤。這里M遠(yuǎn)遠(yuǎn)小于N的維度官地，所以實(shí)現(xiàn)了壓縮。壓縮感知要實(shí)現(xiàn)的是上圖的下半部分烙懦，直接從信號(hào)中采樣M個(gè)成分驱入，然后回復(fù)出近似與原始信號(hào)x的N向量。如果能夠?qū)崿F(xiàn)這一過(guò)程的話氯析，中間的壓縮過(guò)程就沒(méi)有了亏较，將會(huì)極大的方便數(shù)據(jù)存儲(chǔ)和傳輸。
稀疏性/可壓縮性
很多的信號(hào)都是具有稀疏性或者可壓縮性掩缓。Baraniuk的Lecture Note中對(duì)此做了定義雪情。如果存在一個(gè)Nx1維度的信號(hào)x，那么在N維空間中它必然可以表示：
其中你辣，
ψ_i為該空間中N個(gè)正交基向量中的第i個(gè)巡通，并且為N維尘执，對(duì)應(yīng)的s_i為對(duì)應(yīng)正交基向量的系數(shù)。如果一個(gè)信號(hào)x能夠僅用k個(gè)基向量的線性組合來(lái)表示宴凉，那么這個(gè)信號(hào)x被稱為k-sparse誊锭，也就是剩余的N-K個(gè)系數(shù)都為零。可壓縮性和k-sparse的定義稍有不同，如果N個(gè)系數(shù)中眨攘，僅僅只有幾個(gè)數(shù)字非常大郊楣，而其他的數(shù)字則非常小，那么信號(hào)x就具有壓縮性莹捡。
圖像來(lái)源
上圖中一張圖像具有N個(gè)像素，雖然在時(shí)域這N個(gè)像素基本上都不等于零，但是如果對(duì)其進(jìn)行小波變化的話罐盔，基本上除了前面的系數(shù)以為，其他的系數(shù)都很小救崔，對(duì)于聲音信號(hào)等其他信號(hào)惶看，也具有這一特征。也就是說(shuō)六孵，總可以把信號(hào)弄到某個(gè)域中纬黎，讓它在這個(gè)域具有稀疏性。
- 對(duì)于一個(gè)稀疏信號(hào)所在空間劫窒，其存在N個(gè)正交基向量本今，由于其只有K個(gè)位置存在系數(shù)，而其他位置系數(shù)為零主巍，如果把這些點(diǎn)畫(huà)在N維空間中冠息，這些點(diǎn)是多個(gè)K維子空間的集合。比如在三維空間中孕索，有一個(gè)維度的系數(shù)為零逛艰，那么結(jié)果就是幾個(gè)面的組合，如下圖的第一部分搞旭。
- 對(duì)于一個(gè)可壓縮信號(hào)而言散怖，則將構(gòu)建一個(gè)類似于充滿’毛刺‘的子空間，如下圖第二部分所示肄渗。
- 對(duì)于具有稀疏性和可壓縮性的信號(hào)镇眷，他們?cè)诳臻g中的所在位置都是非常靠近空間坐標(biāo)系的恳啥，這是他們共同的特點(diǎn)偏灿。
壓縮感知
假如已知一個(gè)信號(hào)x具有稀疏性，如下圖所示钝的，只有三個(gè)點(diǎn)非零（用綠翁垂、紅铆遭、藍(lán)三個(gè)點(diǎn)表示），如果使用一個(gè)MxN的測(cè)量矩陣（圖中紅色區(qū)域）對(duì)其進(jìn)行測(cè)量沿猜，那么每次的測(cè)量就是x和矩陣某一行向量的內(nèi)積枚荣，得到的測(cè)量結(jié)果就是具有M個(gè)維度的數(shù)據(jù)y。如果M很小的話啼肩，那么就實(shí)現(xiàn)了dimension reduction橄妆。同時(shí)，我們希望M跟K非常的接近祈坠，這樣我們才能夠盡量的保留信號(hào)中的信息害碾。
圖像來(lái)源
一個(gè)很明顯的問(wèn)題是，如果我們的MxN矩陣不是滿秩（full rank）的赦拘，也就是說(shuō)各個(gè)行向量不是獨(dú)立的那么慌随，我們必然會(huì)丟失一些信息。舉個(gè)極端的例子躺同，如果每個(gè)行的行向量是一樣的阁猜，那么得到的y將是M個(gè)相同的值。這顯然是無(wú)法讓我們復(fù)原原始信號(hào)的蹋艺。但是剃袍，如果我們的信號(hào)x是k-sparse的,那么其實(shí)我們只關(guān)心那個(gè)k個(gè)值，對(duì)應(yīng)到矩陣乘法中捎谨，就是我們MxN這個(gè)測(cè)量矩陣中的某K列民效。只要這K列的數(shù)據(jù)是滿秩的就OK了。也就是如果我們能夠找到這樣一個(gè)矩陣涛救，其中任意MxK的子矩陣是full rank研铆，我們就可以用這樣的觀測(cè)矩陣進(jìn)行測(cè)量，并且又可能可以恢復(fù)數(shù)據(jù)州叠。非常幸運(yùn)的是，一個(gè)隨機(jī)的矩陣就有很大的概率具有這樣的性質(zhì)凶赁。所以咧栗，我們可以使用下面的方法來(lái)對(duì)數(shù)據(jù)進(jìn)行隨機(jī)采樣。其中采樣矩陣為MxN維度隨機(jī)矩陣虱肄。采樣得到的每一個(gè)值都是信號(hào)x和隨機(jī)矩陣一個(gè)行向量的內(nèi)積致板。
信號(hào)恢復(fù)
通過(guò)上面的操作，我們獲得了y咏窿，并且也知道自己產(chǎn)生的MxN隨機(jī)矩陣是啥樣斟或，剩下的就是重建我們的信號(hào)x。
因?yàn)榭梢詫?em>x看做在N維空間中集嵌，而這個(gè)MxN矩陣則可以看做一個(gè)在這個(gè)N維空間中的N-M維度的超平面萝挤，這個(gè)超平面穿過(guò)了真正的x御毅，當(dāng)然也穿越了一些別的x，我們要做的就是利用某種標(biāo)準(zhǔn)來(lái)找到真正的x怜珍。大牛們證明L2作為標(biāo)準(zhǔn)是不行滴端蛆，L0雖然是正確的標(biāo)準(zhǔn)，但是很難求解酥泛。如果將L1作為標(biāo)準(zhǔn)也能找到正確的解今豆，但是求解就簡(jiǎn)單很多。
很多的算法能夠支持L1的求解柔袁，但是我們不用關(guān)心這些算法的處理過(guò)程呆躲，只需要對(duì)其進(jìn)行使用即可。
處理在不稀疏的信號(hào)
前面我們提及的幾個(gè)例子中捶索，信號(hào)x只有幾個(gè)地方有值插掂，其余位置的值皆為零。但是如果我們的x是一張圖像的話情组，大部分的位置都將不是零燥筷，這種情況我們要如何處理呢？圖像雖然在時(shí)空域不是稀疏的院崇，但是它在DCT肆氓、小波等變換后的域是稀疏的，所以我可以在這個(gè)域里面處理它底瓣。如下圖所示谢揪，雖然x信號(hào)不稀疏，但是它可以表示為一個(gè)稀疏信號(hào)a和變換矩陣ψ的結(jié)果捐凭。如果將變換矩陣ψ和原來(lái)的矩陣結(jié)合為一個(gè)新的矩陣拨扶，那么就相當(dāng)于在稀疏信號(hào)a上采樣了。
Python實(shí)現(xiàn)壓縮感知
1D信號(hào)
在有了以上的基礎(chǔ)知識(shí)以后茁肠，我們可以來(lái)做一些實(shí)驗(yàn)了患民。主要參考資料來(lái)源Imaging via Compressive Sampling、Compressive Sensing in Python 垦梆。
首先匹颤，我們產(chǎn)生兩個(gè)sin信號(hào)組合而成的信號(hào)
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import scipy.optimize as spopt
import scipy.fftpack as spfft
import scipy.ndimage as spimg
import cvxpy as cvx
n = 5000
t = np.linspace(0, 0.125, n)
y = np.sin(1394 * np.pi * t) + np.sin(3266 * np.pi * t)
yt = spfft.dct(y, norm='ortho')
plt.subplot(221)
plt.plot(t,y,linewidth=0.35)
plt.title('Original signal')
plt.subplot(222)
plt.plot(yt,linewidth=0.5)
plt.title('Frequency Component')
plt.subplot(223)
plt.plot(t[:1000],y[:1000],linewidth=0.45)
plt.title('Original signal, first 1000 points ')
plt.subplot(224)
plt.plot(yt[:500],linewidth=0.5)
plt.title('Original signal, first 1000 points ')
plt.tight_layout()
plt.show()
上面代碼生成的圖像如下：
上面代碼中，出現(xiàn)的下面這一行的目的是對(duì)一維數(shù)據(jù)進(jìn)行DCT變換：
yt = spfft.dct(y, norm='ortho')
可以看到托猩，在DCT域中印蓖，我們的信號(hào)其實(shí)是非常的稀疏的。根據(jù)壓縮感知理論京腥，對(duì)這樣一個(gè)稀疏的信號(hào)赦肃，其實(shí)我們采用隨機(jī)采樣的方式就可以恢復(fù)出原始的信號(hào)。但是，很顯然的一點(diǎn)是他宛，我們無(wú)法直接在DCT域進(jìn)行采樣船侧，只能在時(shí)空域?qū)π盘?hào)進(jìn)行采樣。所以堕汞， 我們對(duì)時(shí)空域的5000個(gè)信號(hào)點(diǎn)進(jìn)行隨機(jī)采樣：
m = 500 # 10% sample
ri = np.random.choice(n, m, replace=False) # random sample of indices
ri.sort() # sorting not strictly necessary, but convenient for plotting
t2 = t[ri]
y2 = y[ri]
plt.plot(t,y,linewidth=0.35)
plt.scatter(t2,y2,marker='o',color='r',s=5)
plt.show()
從上圖中可以看到勺爱，500個(gè)點(diǎn)隨機(jī)分布，雜亂無(wú)章讯检，如何從中提取出原始的信號(hào)呢琐鲁？從
處理在不稀疏的信號(hào)這一節(jié)可以看出，我們可以將這5000個(gè)時(shí)空域的數(shù)據(jù)看做經(jīng)過(guò)反離散余弦變換（iDCT）變換得到的人灼，而"真正的原始信號(hào)"則是‘余弦域’的∥Ф危現(xiàn)在我們要找到一個(gè)變換矩陣，這個(gè)矩陣將余弦域信號(hào)轉(zhuǎn)化為時(shí)空域投放。對(duì)一個(gè)對(duì)角矩陣進(jìn)行反余弦變換奈泪，就能得到這樣一個(gè)矩陣，由于我們對(duì)信號(hào)進(jìn)行了隨機(jī)采樣灸芳，所以將隨機(jī)采樣的位置也考慮其中涝桅。最后得到的A矩陣就是我們整體的采樣矩陣，包含了變換和隨機(jī)采樣的信息烙样。
A = spfft.idct(np.identity(n), norm='ortho', axis=0)
A = A[ri]
上面的第一行冯遂，首先生成了這樣一個(gè)矩陣，然后第二行按照隨機(jī)采樣的方式取出相應(yīng)的行谒获。
到這里蛤肌，我們做好了所有的準(zhǔn)備工作：
- 得到了采樣后的數(shù)據(jù)點(diǎn)
- 得到了采樣矩陣
剩下的就是在L1條件下，求解出我們的原始信號(hào)批狱。然后在對(duì)得到的信號(hào)進(jìn)行反余弦變換裸准，就得到了時(shí)空域的信號(hào)了。
vx = cvx.Variable(n)
objective = cvx.Minimize(cvx.norm(vx, 1))
constraints = [A*vx == y2]
prob = cvx.Problem(objective, constraints)
result = prob.solve(verbose=True)
plt.plot(sig,linewidth=0.35)
plt.show()
2D信號(hào)
在重建了一個(gè)1D信號(hào)以后赔硫，我們可以對(duì)2D信號(hào)進(jìn)行壓縮采樣并重構(gòu)炒俱。很顯然，一般的圖像信號(hào)在時(shí)空域也是非稀疏的爪膊，我們?nèi)匀恍枰O(shè)計(jì)一個(gè)轉(zhuǎn)換矩陣向胡，將其轉(zhuǎn)化到一個(gè)稀疏域下，這里惊完，我們繼續(xù)使用DCT變換。根據(jù)2D-Signal Transforms and The Separability Property中的第7頁(yè)的內(nèi)容处硬，可以看到小槐，對(duì)于一個(gè)可分的圖像變換，可以將其看做兩個(gè)一維變換分別作用于圖像矩陣的行和列。這一點(diǎn)類似高斯濾波的可分性凿跳。而行變換和列變換的克羅內(nèi)克積則等效于這個(gè)2D變換件豌。也就是說(shuō)，如果將圖像展開(kāi)為一維數(shù)據(jù)以后控嗜，乘以這個(gè)克羅內(nèi)克積得到的矩陣茧彤，就相當(dāng)于對(duì)圖像做了2D iDCT變換：
下面我們來(lái)進(jìn)行對(duì)圖像隨機(jī)采樣，并且復(fù)原的工作疆栏，首先曾掂，我們讀取lena圖，然后隨機(jī)取出里面50%的像素
Xorig = spimg.imread('lena_gray.gif',flatten=True, mode='L')
X = spimg.zoom(Xorig,0.25)
ny,nx=X.shape
k = round(nx*ny*0.5)
ri = np.random.choice(nx*ny,np.int(k), replace=False)
plt.subplot(121)
plt.imshow(Xorig,cmap='gray')
plt.axis('off')
plt.title('Original')
plt.subplot(122)
X1d = X.flatten()
b = X1d[ri]
Xsam = np.zeros(np.shape(X.flatten()))
Xsam[ri] = b;
Xsam = np.reshape(Xsam,[ny,nx])
plt.imshow(Xsam,cmap='gray')
plt.title('50% Sample')
plt.axis('off')
plt.show()
從上面的結(jié)果可以看到壁顶，通過(guò)對(duì)Lena圖像進(jìn)行50%隨機(jī)的采樣珠洗，最后得到的圖像已經(jīng)模糊不清了，但是我們?nèi)搜廴匀皇强梢钥闯鰞?nèi)容是Lena圖像若专。
A = np.kron(
spfft.idct(np.identity(nx), norm='ortho', axis=0),
spfft.idct(np.identity(ny), norm='ortho', axis=0)
)
