學習壓縮感知比較好的文章鏈接收藏

一、首推：壓縮感知“Hello World”代碼初步學習

本文真的對我理解壓感起到很大的幫助，特別在OMP（正交匹配追蹤重構算法）算法的講解中，尤為精彩，我從來沒有見過哪個博主能如此認真的去一行一行地解釋，解讀代碼，我對此十分感謝。

個人認為，能有超過70%注釋的代碼才是負責人的代碼，當然也是好的代碼。本文做到的不僅如此，簡直優秀。為此，我一行一行的抄寫了里面的代碼。

精彩節選：

1. 原始信號x是什么？我采集的是原始信號x還是y = Ax得到的y？

解答：

記原始信號為x，我們在sensor方得到的原始信號就是n*1的信號x，而在receiver方采集到的信號是y。針對y=Ax做變換時，A（m*n ）是一個隨機矩陣（真的很隨機，不用任何正交啊什么的限定）。通過由隨機矩陣變換內積得到y，我們的目標是從y中恢復x。由于A是m*n(m

（我認為藍色部分是精髓）

2. O M P 重構算法：

% 1-D信號壓縮傳感的實現(正交匹配追蹤法Orthogonal Matching Pursuit)

% 測量數M>=K*log(N/K),K是稀疏度,N信號長度,可以近乎完全重構

% 編程人--香港大學電子工程系 沙威 Email: wsha@eee.hku.hk

% 編程時間：2008年11月18日

% 文檔下載: http://www.eee.hku.hk/~wsha/Freecode/freecode.htm

% 參考文獻：Joel A. Tropp and Anna C. Gilbert

% Signal Recovery From Random Measurements Via Orthogonal Matching

% Pursuit，IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,

% DECEMBER 2007.

clc;clear

%% 1. 時域測試信號生成

K=7; % 稀疏度(做FFT可以看出來)

N=256; % 信號長度

M=64; % 測量數(M>=K*log(N/K),至少40,但有出錯的概率)

f1=50; % 信號頻率1

f2=100; % 信號頻率2

f3=200; % 信號頻率3

f4=400; % 信號頻率4

fs=800; % 采樣頻率

ts=1/fs; % 采樣間隔

Ts=1:N; % 采樣序列

x=0.3*cos(2*pi*f1*Ts*ts)+0.6*cos(2*pi*f2*Ts*ts)+0.1*cos(2*pi*f3*Ts*ts)+0.9*cos(2*pi*f4*Ts*ts); % 完整信號,由4個信號疊加而來

%% 2. 時域信號壓縮傳感

Phi=randn(M,N); % 測量矩陣(高斯分布白噪聲)64*256的扁矩陣，Phi也就是文中說的D矩陣

s=Phi*x.'; % 獲得線性測量 ，s相當于文中的y矩陣

%% 3. 正交匹配追蹤法重構信號(本質上是L_1范數最優化問題)

%匹配追蹤：找到一個其標記看上去與收集到的數據相關的小波；在數據中去除這個標記的所有印跡；不斷重復直到我們能用小波標記“解釋”收集到的所有數據。

m=2*K; % 算法迭代次數(m>=K)，設x是K-sparse的

Psi=fft(eye(N,N))/sqrt(N); % 傅里葉正變換矩陣

T=Phi*Psi'; % 恢復矩陣(測量矩陣*正交反變換矩陣)

hat_y=zeros(1,N); % 待重構的譜域(變換域)向量

Aug_t=[]; % 增量矩陣(初始值為空矩陣)

r_n=s; % 殘差值

for times=1:m; % 迭代次數(有噪聲的情況下,該迭代次數為K)

for col=1:N; % 恢復矩陣的所有列向量

product(col)=abs(T(:,col)'*r_n); % 恢復矩陣的列向量和殘差的投影系數(內積值)

end

[val,pos]=max(product); % 最大投影系數對應的位置，即找到一個其標記看上去與收集到的數據相關的小波

Aug_t=[Aug_t,T(:,pos)]; % 矩陣擴充

T(:,pos)=zeros(M,1); % 選中的列置零（實質上應該去掉，為了簡單我把它置零），在數據中去除這個標記的所有印跡

aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s; % 最小二乘,使殘差最小

r_n=s-Aug_t*aug_y; % 殘差

pos_array(times)=pos; % 紀錄最大投影系數的位置

end

hat_y(pos_array)=aug_y; % 重構的譜域向量

hat_x=real(Psi'*hat_y.'); % 做逆傅里葉變換重構得到時域信號

%% 4. 恢復信號和原始信號對比

figure(1);

hold on;

plot(hat_x,'k.-') % 重建信號

plot(x,'r') % 原始信號

legend('Recovery','Original')

norm(hat_x.'-x)/norm(x)

3.對上述代碼的解讀：

CS的前提是信號的稀疏性，這包括信號本身在時域上是稀疏的或者信號經過一定的變換在相應的變換域（包括頻域、小波域等）上是稀疏的。通過y=Phi*x得到測量信號（y是M維，Phi是M*N維測量矩陣，x為N維原始信號），這樣得到的測量信號y就可以傳輸或者保存了。但是CS的關鍵問題是通過一定的算法由y恢復原始信號x，恢復方法是先由y得到原始信號x在變換域上的估計值hat_y，得到hat_y之后經過反變換自然就得到了原始信號的估計值hat_x。本文中講到的OMP方法實際上是解決如何由y得到hat_y。

最精彩部分：

好了，有了OMP算法，開始對應解釋代碼：

for col=1:N; ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?% ?恢復矩陣的所有列向量

product(col)=abs(T(:,col)'*r_n); ? ? ? ? ?% ?恢復矩陣的列向量和殘差的投影系數(內積值)

end

這個循環是讓矩陣T的每一列與殘差求內各，T一共有N列，這里得到N個內積值存在product里面。內積值最大的即為相關性最強T(:,col)為M*1列向量，r_n初如化為s，是M*1列向量，這里讓T(:,col)轉置后再與r_n相乘，即一個1*M的行向量與一個M*1的列向量相乘，根據矩陣運算規則結果為一個數（即1*1的矩陣）。

[val,pos]=max(product); 這句話的關鍵是得到pos，即得到T中的哪一列與殘差r_n的內積值最大，也就是哪一列與殘差r_n相關性最強。此即 英文步驟中的第二步 。

Aug_t=[Aug_t,T(:,pos)]; 此即 英文步驟中的第三步 ，將剛剛得到的與殘差r_n內積值最大的列存到Aug_t中，這個矩陣隨著循環次數（迭代次數）的變換而變化，是M*times的矩陣。

T(:,pos)=zeros(M,1); 這一句是為了下一次迭代做準備的，這次找到了與殘差最相關的列，將殘差更新后，下次再找與殘差僅次于這一列的T的另外一列；

aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s; 這一句即 英文步驟中的第四步 ，這句加上后面一句也是困擾了我好久兩句代碼，所以得說說：

首先我們針對的是s=T*hat_y，現在是已知s要求hat_y，現在假如說矩陣T為N*N方陣且滿秩（即N個未知數，N個獨立的方程），那么很容易知道hat_y=T^-1 * s，其中T^-1表示矩陣T的逆矩陣。但是現在T是一個M*N的扁矩陣，矩陣T沒有常規意義上的逆矩陣，這里就有“廣義逆”的概念（詳情參見國內矩陣分析教材），hat_y的解可能是不存在的，我們這里要求的是 最小二乘解 aug_y，最小二乘解aug_y將使s-T*aug_y這個列向量2范數最小。

對于用矩陣形式表達的線性方程組：

它的最小二乘解為：

其中

即為矩陣G的最小二乘廣義逆（廣義逆的一種）。

有了這些知識背景后代碼就容易理解了，在第三步中，得到矩陣T中的與殘差r_n最相關的列組成的矩陣Aug_t，而第四步實際上就是在求方程組Aug_t*Aug_y=s的最小二乘解。

r_n=s-Aug_t*aug_y;這一句就是用求得的最小二乘解更新殘差r_n，在下一次迭代中使用。注意最小二乘解的含義，它并不是使Aug_t*Aug_y=s成立，而只是讓s-Aug_t*aug_y的2范數最小，而r_n就是最小的值。此即 英文步驟中的第五步 ，兩個式子合在一起寫了。

pos_array(times)=pos; 把與T中與殘差最相關的列號記下來，恢復時使用。

到此，主要的for循環就說完了。

hat_y(pos_array)=aug_y; 最后一次迭代得到的最小二乘解aug_y即為恢復的值，位置分別對應于迭代中每一次與殘差r_s最相關的矩陣T的列號。hat_y(pos_array)大小是和pos_array大小一樣的，并且 hat_y(pos_array)的第k個元素就是pos_array(k)。

hat_x=real(Psi'*hat_y.');此即：

，這里用hat_x以與原如信號x區分，x為原信號，hat_x為恢復的信號。代碼中對hat_y取了轉置是因為hat_y應該是個列向量，而在代碼中的前面hat_y=zeros(1,N); 將其命成了行向量，所以這里轉置了一下，沒什么大不了的。

matlab運行一下就發現了 psi*hat_y的結果是實數，但是都帶著+0.0000i 所以要取實部。

————————————————————————————————————————————————————

是不是對這段代碼解讀感動的淚流滿面。

事實上，學習這段代碼的目的并不在此例本身，而是通過對這段代碼的學習，來舉一反三，應用于自己的工程，或者用于理解自己的工程相關的代碼。

為什么會感動，那是因為看到自己的前輩留下來的注釋率不到2%的代碼砸到你的臉上讓你死磕，那種無奈與孤獨可曾了解？看到這種代碼能不感動嗎？這是我等的典范！！！

————————————————————————————————————————————————————

二、姊妹篇，輔助理解上述文章的一篇博文《“壓縮感知” 之 “Hello World”》

起名為：Hello World，我服。這才是配得上Hello world 之名的博文。

三、基于壓縮傳感的匹配追蹤重構算法研究

這是一篇碩士學位論文，對于研究壓感重建算法很有幫助。

主要是中文版的，所以閱讀起來障礙會小一點。

文中給出了O M P重建算法的中文版流程：

這只是本論文中的冰山一角，論文繼續研讀中。

四、寬帶接收機中的非均勻采樣技術研究

這是西安電子科技大學的一篇碩士論文，哈哈，這是我師兄的一篇畢業碩士畢業論文，雖然沒有見過師兄，但覺得師兄應該挺牛的，這是第一個研究這個項目的人，目前此項目由我接手，說來已經經歷第三代實驗室人了。同時也感謝親手將此項目傳于我的師兄，讓我更快的走上了硬件的這條路，還有我的即將離開的研三師兄，你對我的指導也讓我感激不盡。

持續更新中。

下面一篇文章，我即將閱讀，先保存于此：

MP算法和OMP算法及其思想

機器學習

版權聲明：本文內容由網絡用戶投稿，版權歸原作者所有，本站不擁有其著作權，亦不承擔相應法律責任。如果您發現本站中有涉嫌抄襲或描述失實的內容，請聯系我們jiasou666@gmail.com 處理，核實后本網站將在24小時內刪除侵權內容。