一、首推：壓縮感覺“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<n)的，是以原信号f（n*1）信号被壓縮到y（m*1）。
（我認為藍色部分是精髓）
2. O M P 重構算法：
% 1-D信号壓縮傳感的實作(正交比對追蹤法Orthogonal Matching Pursuit)
% 測量數M>=K*log(N/K),K是稀疏度,N信号長度,可以近乎完全重構
% 程式設計人--香港大學電子工程系 沙威 Email: [email protected]
% 程式設計時間：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重建算法的中文版流程：
這隻是本論文中的冰山一角，論文繼續研讀中。
四、寬帶接收機中的非均勻采樣技術研究
這是西安電子科技大學的一篇碩士論文，哈哈，這是我師兄的一篇畢業碩士畢業論文，雖然沒有見過師兄，但覺得師兄應該挺牛的，這是第一個研究這個項目的人，目前此項目由我接手，說來已經經曆第三代實驗室人了。同時也感謝親手将此項目傳于我的師兄，讓我更快的走上了硬體的這條路，還有我的即将離開的研三師兄，你對我的指導也讓我感激不盡。
持續更新中。