正交匹配追蹤(OMP)算法的MATLAB函數代碼並給出單次測試例程代碼
測量數M與重構成功概率關系曲線繪制例程代碼
信號稀疏度K與重構成功概率關系曲線繪制例程代碼
參考來源：http://blog.csdn.net/jbb0523/article/details/45130793
參考文獻：Joel A. Tropp and Anna C. Gilbert. Signal Recovery From Random Measurements Via Orthogonal Matching Pursuit[J]. IEEETransactions on Information Theory, VOL. 53, NO. 12, DECEMBER 2007.
0、符號說明如下
壓縮觀測y=Φx，其中y為觀測所得向量M×1，x為原信號N×1（M<<N）。x一般不是稀疏的，但在某個變換域Ψ是稀疏的，即x=Ψθ，其中θ為K稀疏的，即θ只有K個非零項。此時y=ΦΨθ，令A=ΦΨ，則y=Aθ。
(1)y為觀測所得向量，大小為M×1
(2)x為原信號，大小為N×1
(3)θ為K稀疏的，是信號在x在某變換域的稀疏表示
(4)Φ稱為觀測矩陣、測量矩陣、測量基，大小為M×N
(5)Ψ稱為變換矩陣、變換基、稀疏矩陣、稀疏基、正交基字典矩陣，大小為N×N
(6)A稱為測度矩陣、傳感矩陣、CS信息算子，大小為M×N
上式中，一般有K<<M<<N，后面三個矩陣各個文獻的叫法不一，以后我將Φ稱為測量矩陣、將Ψ稱為稀疏矩陣、將A稱為傳感矩陣。
1、OMP重構算法流程
2、正交匹配追蹤(OMP)MATLAB代碼(CS_OMP.m)
function[theta]=CS_OMP(y,A,t)
%CS_OMP Summary of this function goes here
%Version: 1.0 written by jbb0523 @2015-04-18
% Detailed explanation goes here
% y = Phi * x
% x = Psi * theta
% y = Phi*Psi * theta
% 令 A = Phi*Psi, 則y=A*theta
% 現在已知y和A，求theta
[y_rows,y_columns]=size(y);
ify_rows<y_columns
y=y';%y should be a column vector
end
[M,N]=size(A); %傳感矩陣A為M*N矩陣
theta=zeros(N,1); %用來存儲恢復的theta(列向量)
At=zeros(M,t); %用來迭代過程中存儲A被選擇的列
Pos_theta=zeros(1,t); %用來迭代過程中存儲A被選擇的列序號
r_n=y; %初始化殘差(residual)為y
forii=1:t %迭代t次，t為輸入參數
product=A'*r_n; %傳感矩陣A各列與殘差的內積
[val,pos]=max(abs(product)); %找到最大內積絕對值，即與殘差最相關的列
At(:,ii)=A(:,pos); %存儲這一列
Pos_theta(ii)=pos; %存儲這一列的序號
A(:,pos)=zeros(M,1); %清零A的這一列，其實此行可以不要，因為它與殘差正交
%y=At(:,1:ii)*theta，以下求theta的最小二乘解(Least Square)
theta_ls=(At(:,1:ii)'*At(:,1:ii))^(-1)*At(:,1:ii)'*y; %最小二乘解
%At(:,1:ii)*theta_ls是y在At(:,1:ii)列空間上的正交投影
r_n=y-At(:,1:ii)*theta_ls; %更新殘差
end
theta(Pos_theta)=theta_ls; %恢復出的theta
end
3、OMP單次重構測試代碼(CS_Reconstuction_Test.m)
代碼中，直接構造一個K稀疏的信號，所以稀疏矩陣為單位陣。
%壓縮感知重構算法測試
clear all;close all;clc;
M=64;%觀測值個數
N=256;%信號x的長度
K=10;%信號x的稀疏度
Index_K=randperm(N);
x=zeros(N,1);
x(Index_K(1:K))=5*randn(K,1);%x為K稀疏的，且位置是隨機的
Psi=eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta
Phi=randn(M,N);%測量矩陣為高斯矩陣
A=Phi*Psi;%傳感矩陣
y=Phi*x;%得到觀測向量y
%% 恢復重構信號x
tic
theta=CS_OMP(y,A,K);
x_r=Psi*theta;% x=Psi * theta
toc
%% 繪圖
figure;
plot(x_r,'k.-');%繪出x的恢復信號
hold on;
plot(x,'r');%繪出原信號x
hold off;
legend('Recovery','Original')
fprintf('\n恢復殘差：');
norm(x_r-x)%恢復殘差
代碼解釋：上述代碼是直接構造一個K稀疏的信號，接下來解釋一下代碼中是如何構造該稀疏信號的。
首先介紹一下randperm函數，即randm permutation隨機排列、隨機置換。功能是隨機打亂一個數字序列，其內的參數決定了隨機數的范圍。
Index_K = randperm(N); 指的是將1到256的數進行隨機排列
初始化信號x：x = zeros(N,1);
x(Index_K(1:K)) = 5*randn(K,1)；等式右邊比較好理解，即構造一個K*1的隨機向量，接着解釋等式左邊，括號內Index_K（1:K）指的是選取隨機排列后的數列的前K項，因為我們要構造的信號是K稀疏的，也就是只有K個項為非零元素。則我們要將等式右邊產生的K個非零值隨機的插到信號x的K個位置中，舉個例子，比如經過排列后的Index_K（1:K）=12 56 30 17 5 2 6 98 200 85 ，則等式右邊的K個非零值被放置在x的第12、56……的位置上。
接着解釋最后一行代碼，norm指的是范數的意思，在代碼中求得是重構后的信號與原始信號的差值的一范數，一范數相當於求絕對值，據此求出誤差。
運行結果如下：（信號為隨機生成，所以每次結果均不一樣）
1）圖：
2）Command Windows
Elapsed time is 0.849710 seconds.
恢復殘差：
ans =
5.5020e-015
4、測量數M與重構成功概率關系曲線繪制例程代碼
這段代碼真的是斷斷續續看了好久才明白，理解代碼還是要分塊分塊搞懂
%壓縮感知重構算法測試CS_Reconstuction_MtoPercentage.m
% 繪制參考文獻中的Fig.1
% 參考文獻：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.
% Elapsed time is 1171.606254 seconds.(@20150418night)
clear all;close all;clc;
%% 參數配置初始化
CNT=1000;%對於每組(K,M,N)，重復迭代次數
N=256;%信號x的長度
Psi=eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta
K_set=[4,12,20,28,36];%信號x的稀疏度集合
Percentage=zeros(length(K_set),N);%存儲恢復成功概率
%% 主循環，遍歷每組(K,M,N)
tic
forkk=1:5
K=K_set(kk);%本次稀疏度
M_set=K:5:N;%M沒必要全部遍歷，每隔5測試一個就可以了
PercentageK=zeros(1,length(M_set));%存儲此稀疏度K下不同M的恢復成功概率
formm=1:length(M_set)
M=M_set(mm);%本次觀測值個數
P=0;
forcnt=1:CNT%每個觀測值個數均運行CNT次
Index_K=randperm(N);
x=zeros(N,1);
x(Index_K(1:K))=5*randn(K,1);%x為K稀疏的，且位置是隨機的
Phi=randn(M,N);%測量矩陣為高斯矩陣
A=Phi*Psi;%傳感矩陣
y=Phi*x;%得到觀測向量y
theta=CS_OMP(y,A,K);%恢復重構信號theta
x_r=Psi*theta;% x=Psi * theta
ifnorm(x_r-x)<1e-6%如果殘差小於1e-6則認為恢復成功
P=P+1;
end
end
PercentageK(mm)=P/CNT*100;%計算恢復概率
end
Percentage(kk,1:length(M_set))=PercentageK;
end
toc
save MtoPercentage1000%運行一次不容易，把變量全部存儲下來
%% 繪圖
S=['-ks';'-ko';'-kd';'-kv';'-k*'];
figure;
forkk=1:length(K_set)
K=K_set(kk);
M_set=K:5:N;
L_Mset=length(M_set);
plot(M_set,Percentage(kk,1:L_Mset),S(kk,:));%繪出x的恢復信號
hold on;
end
hold off;
xlim([0256]);
legend('K=4','K=12','K=20','K=28','K=36');
xlabel('Number of measurements(M)');
ylabel('Percentage recovered');
title('Percentage of input signals recovered correctly(N=256)(Gaussian)');
這段代碼指的是固定稀疏度的情況下，研究測量次數與重構概率的關系，突然間不知道測量次數指什么。觀測矩陣大小為M*N，測量次數也就是指的這個M的大小，也就是我們壓縮采樣的采樣點數。
8-14行代碼都是初始化，但是第14行代碼，Percentage = zeros(length(K_set),N);為什么維度要這樣設置呢？我們要得出的圖形是以測量次數M為橫坐標，重構概率為縱坐標的，測量次數最大為數據的長度，也就是N，因為我們在仿真中對不同稀疏度的情況進行了仿真，共仿真5種不同稀疏度的情況，所以行數為5，即length（K_set）
接着在第17行進入了主循環，第19行M_set = K:5:N；沒必要全部遍歷，所以每隔5個對該點的值進行測試，但為什么要從K開始呢？K指的是信號的稀疏度，就是信號x最多的非零元素，所以我們進行觀測的時候最少要觀測到所有非零元素，所以從K開始。執行完這行代碼之后生成一個測量次數的行向量，注意不同稀疏度下的測量次數集合是不同的。
選擇了此次測試的稀疏度后，第21行代碼開始對該稀疏度下的測量次數與重構精度的關系進行了測試。依次 選擇測量次數集合M_set中的測量次數，第23行初始化P=0，后面如果殘差小於某一個值時，即重構成功時，P+1。每個觀測值重復1000次操作。
第25到32行是生成稀疏信號並進行OMP重構，得到重構后的信號。
第37行代碼，重復試驗1000次后，記錄下當前測量次數下的恢復概率，P指的是重構成功的個數，除以1000次試驗次數再乘上100即得到重構的概率。
接着進行下一個觀測次數的循環。
M_set集合中的測量次數全部執行完畢后，執行第39行代碼：Percentage(kk,1:length(M_set)) = PercentageK，將此次稀疏度下測得的重構概率保存到Percentage中，Percentage的行數是稀疏度的個數，列數是測量次數的個數。
第44行代碼開始是繪圖，根據稀疏度先得到測量次數的集合，然后以測量次數M為橫軸，重構概率為縱軸繪制圖形。
本程序運行結果：
文獻中的Fig.1：
5、信號稀疏度K與重構成功概率關系曲線繪制例程代碼
代碼與4中的類似
%壓縮感知重構算法測試CS_Reconstuction_KtoPercentage.m % 繪制參考文獻中的Fig.2 % 參考文獻：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. % Elapsed time is 1448.966882 seconds.(@20150418night) clear all;close all;clc; %% 參數配置初始化 CNT=1000;%對於每組(K,M,N)，重復迭代次數 N=256;%信號x的長度 Psi=eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta M_set=[52,100,148,196,244];%測量值集合 Percentage=zeros(length(M_set),N);%存儲恢復成功概率 %% 主循環，遍歷每組(K,M,N) tic formm=1:length(M_set) M=M_set(mm);%本次測量值個數 K_set=1:5:ceil(M/2);%信號x的稀疏度K沒必要全部遍歷，每隔5測試一個就可以了 PercentageM=zeros(1,length(K_set));%存儲此測量值M下不同K的恢復成功概率 forkk=1:length(K_set) K=K_set(kk);%本次信號x的稀疏度K P=0; forcnt=1:CNT%每個觀測值個數均運行CNT次 Index_K=randperm(N); x=zeros(N,1); x(Index_K(1:K))=5*randn(K,1);%x為K稀疏的，且位置是隨機的 Phi=randn(M,N);%測量矩陣為高斯矩陣 A=Phi*Psi;%傳感矩陣 y=Phi*x;%得到觀測向量y theta=CS_OMP(y,A,K);%恢復重構信號theta x_r=Psi*theta;% x=Psi * theta ifnorm(x_r-x)<1e-6%如果殘差小於1e-6則認為恢復成功 P=P+1; end end PercentageM(kk)=P/CNT*100;%計算恢復概率 end Percentage(mm,1:length(K_set))=PercentageM; end toc save KtoPercentage1000test%運行一次不容易，把變量全部存儲下來 %% 繪圖 S=['-ks';'-ko';'-kd';'-kv';'-k*']; figure; formm=1:length(M_set) M=M_set(mm); K_set=1:5:ceil(M/2); L_Kset=length(K_set); plot(K_set,Percentage(mm,1:L_Kset),S(mm,:));%繪出x的恢復信號 hold on; end hold off; xlim([0125]); legend('M=52','M=100','M=148','M=196','M=244'); xlabel('Sparsity level(K)'); ylabel('Percentage recovered'); title('Percentage of input signals recovered correctly(N=256)(Gaussian)');
本程序運行結果：
文獻中的Fig.2：
