壓縮感知重構(gòu)算法之廣義正交匹配追蹤

主要內(nèi)容：

gOMP的算法流程gOMP的MATLAB實(shí)現(xiàn)一維信號(hào)的實(shí)驗(yàn)與結(jié)果稀疏度K與重構(gòu)成功概率關(guān)系的實(shí)驗(yàn)與結(jié)果。

一、gOMP的算法流程

廣義正交匹配追蹤(Generalized OMP, gOMP)算法可以看作為OMP算法的一種推廣。OMP每次只選擇與殘差相關(guān)最大的一個(gè)，而gOMP則是簡(jiǎn)單地選擇最大的S個(gè)。之所以這里表述為"簡(jiǎn)單地選擇"是相比于ROMP之類(lèi)算法的，不進(jìn)行任何其它處理，只是選擇最大的S個(gè)而已。

gOMP的算法流程：

二、gOMP的MATLAB實(shí)現(xiàn)(CS_gOMP.m)

function?[?theta?]?=?CS_gOMP(?y,A,K,S?)
%???CS_gOMP
%???Detailed?explanation?goes?here
%???y?=?Phi?*?x
%???x?=?Psi?*?theta
%????y?=?Phi*Psi?*?theta
%???令?A?=?Phi*Psi,?則y=A*theta
%???現(xiàn)在已知y和A，求theta
%???Reference:?Jian?Wang,?Seokbeop?Kwon,?Byonghyo?Shim.??Generalized?
%???orthogonal?matching?pursuit,?IEEE?Transactions?on?Signal?Processing,?
%???vol.?60,?no.?12,?pp.?6202-6216,?Dec.?2012.?
%???Available?at:?http://islab.snu.ac.kr/paper/tsp_gOMP.pdf
????if?nargin?<?4
????????S?=?round(max(K/4,?1));
????end
????[y_rows,y_columns]?=?size(y);
????if?y_rowsM
????????????if?ii?==?1
????????????????theta_ls?=?0;
????????????end
????????????break;
????????end
????????At?=?A(:,Sk);%將A的這幾列組成矩陣At
????????%y=At*theta，以下求theta的最小二乘解(Least?Square)
????????theta_ls?=?(At'*At)^(-1)*At'*y;%最小二乘解
????????%At*theta_ls是y在At)列空間上的正交投影
????????r_n?=?y?-?At*theta_ls;%更新殘差
????????Pos_theta?=?Sk;
????????if?norm(r_n)<1e-6
????????????break;%quit?the?iteration
????????end
????end
????theta(Pos_theta)=theta_ls;%恢復(fù)出的theta
end

三、一維信號(hào)的實(shí)驗(yàn)與結(jié)果

%壓縮感知重構(gòu)算法測(cè)試
clear?all;close?all;clc;
M?=?128;%觀測(cè)值個(gè)數(shù)
N?=?256;%信號(hào)x的長(zhǎng)度
K?=?30;%信號(hào)x的稀疏度
Index_K?=?randperm(N);
x?=?zeros(N,1);
x(Index_K(1:K))?=?5*randn(K,1);%x為K稀疏的，且位置是隨機(jī)的
Psi?=?eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta
Phi?=?randn(M,N)/sqrt(M);%測(cè)量矩陣為高斯矩陣
A?=?Phi?*?Psi;%傳感矩陣
y?=?Phi?*?x;%得到觀測(cè)向量y

%%?恢復(fù)重構(gòu)信號(hào)x
tic
theta?=?CS_gOMP(?y,A,K);
x_r?=?Psi?*?theta;%?x=Psi?*?theta
toc

%%?繪圖
figure;
plot(x_r,'k.-');%繪出x的恢復(fù)信號(hào)
hold?on;
plot(x,'r');%繪出原信號(hào)x
hold?off;
legend('Recovery','Original')
fprintf('n恢復(fù)殘差：');
norm(x_r-x)%恢復(fù)殘差

四、稀疏數(shù)K與重構(gòu)成功概率關(guān)系的實(shí)驗(yàn)與結(jié)果

%???壓縮感知重構(gòu)算法測(cè)試CS_Reconstuction_KtoPercentagegOMP.m
%???Reference:?Jian?Wang,?Seokbeop?Kwon,?Byonghyo?Shim.??Generalized?
%???orthogonal?matching?pursuit,?IEEE?Transactions?on?Signal?Processing,?
%???vol.?60,?no.?12,?pp.?6202-6216,?Dec.?2012.?
%???Available?at:?http://islab.snu.ac.kr/paper/tsp_gOMP.pdf

clear?all;close?all;clc;
addpath(genpath('../../OMP/'))

%%?參數(shù)配置初始化
CNT?=?1000;?%對(duì)于每組(K,M,N)，重復(fù)迭代次數(shù)
N?=?256;?%信號(hào)x的長(zhǎng)度
Psi?=?eye(N);?%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta
M_set?=?[128];?%測(cè)量值集合
KIND?=?['OMP??????';'ROMP?????';'StOMP????';'SP???????';'CoSaMP???';...
????'gOMP(s=3)';'gOMP(s=6)';'gOMP(s=9)'];
Percentage?=?zeros(N,length(M_set),size(KIND,1));?%存儲(chǔ)恢復(fù)成功概率

%%?主循環(huán)，遍歷每組(K,M,N)
tic
for?mm?=?1:length(M_set)
????M?=?M_set(mm);?%本次測(cè)量值個(gè)數(shù)
????K_set?=?5:5:70;?%信號(hào)x的稀疏度K沒(méi)必要全部遍歷，每隔5測(cè)試一個(gè)就可以了
????%存儲(chǔ)此測(cè)量值M下不同K的恢復(fù)成功概率
????PercentageM?=?zeros(size(KIND,1),length(K_set));
????for?kk?=?1:length(K_set)
???????K?=?K_set(kk);?%本次信號(hào)x的稀疏度K
???????P?=?zeros(1,size(KIND,1));
???????fprintf('M=%d,K=%dn',M,K);
???????for?cnt?=?1:CNT??%每個(gè)觀測(cè)值個(gè)數(shù)均運(yùn)行CNT次
????????????Index_K?=?randperm(N);
????????????x?=?zeros(N,1);
????????????x(Index_K(1:K))?=?5*randn(K,1);?%x為K稀疏的，且位置是隨機(jī)的????????????????
????????????Phi?=?randn(M,N)/sqrt(M);?%測(cè)量矩陣為高斯矩陣
????????????A?=?Phi?*?Psi;?%傳感矩陣
????????????y?=?Phi?*?x;?%得到觀測(cè)向量y
????????????%(1)OMP
????????????theta?=?CS_OMP(y,A,K);?%恢復(fù)重構(gòu)信號(hào)theta
????????????x_r?=?Psi?*?theta;?%?x=Psi?*?theta
????????????if?norm(x_r-x)<1e-6?%如果殘差小于1e-6則認(rèn)為恢復(fù)成功
????????????????P(1)?=?P(1)?+?1;
????????????end
????????????%(2)ROMP
????????????theta?=?CS_ROMP(y,A,K);?%恢復(fù)重構(gòu)信號(hào)theta
????????????x_r?=?Psi?*?theta;?%?x=Psi?*?theta
????????????if?norm(x_r-x)<1e-6?%如果殘差小于1e-6則認(rèn)為恢復(fù)成功
????????????????P(2)?=?P(2)?+?1;
????????????end
????????????%(3)StOMP
????????????theta?=?CS_StOMP(y,A);?%恢復(fù)重構(gòu)信號(hào)theta
????????????x_r?=?Psi?*?theta;?%?x=Psi?*?theta
????????????if?norm(x_r-x)<1e-6?%如果殘差小于1e-6則認(rèn)為恢復(fù)成功
????????????????P(3)?=?P(3)?+?1;
????????????end
????????????%(4)SP
????????????theta?=?CS_SP(y,A,K);?%恢復(fù)重構(gòu)信號(hào)theta
????????????x_r?=?Psi?*?theta;?%?x=Psi?*?theta
????????????if?norm(x_r-x)<1e-6?%如果殘差小于1e-6則認(rèn)為恢復(fù)成功
????????????????P(4)?=?P(4)?+?1;
????????????end
????????????%(5)CoSaMP
????????????theta?=?CS_CoSaMP(y,A,K);?%恢復(fù)重構(gòu)信號(hào)theta
????????????x_r?=?Psi?*?theta;?%?x=Psi?*?theta
????????????if?norm(x_r-x)<1e-6?%如果殘差小于1e-6則認(rèn)為恢復(fù)成功
????????????????P(5)?=?P(5)?+?1;
????????????end
????????????%(6)gOMP,S=3
????????????theta?=?CS_gOMP(y,A,K,3);?%恢復(fù)重構(gòu)信號(hào)theta
????????????x_r?=?Psi?*?theta;?%?x=Psi?*?theta
????????????if?norm(x_r-x)<1e-6?%如果殘差小于1e-6則認(rèn)為恢復(fù)成功
????????????????P(6)?=?P(6)?+?1;
????????????end
????????????%(7)gOMP,S=6
????????????theta?=?CS_gOMP(y,A,K,6);?%恢復(fù)重構(gòu)信號(hào)theta
????????????x_r?=?Psi?*?theta;?%?x=Psi?*?theta
????????????if?norm(x_r-x)<1e-6?%如果殘差小于1e-6則認(rèn)為恢復(fù)成功
????????????????P(7)?=?P(7)?+?1;
????????????end
????????????%(8)gOMP,S=9
????????????theta?=?CS_gOMP(y,A,K,9);?%恢復(fù)重構(gòu)信號(hào)theta
????????????x_r?=?Psi?*?theta;?%?x=Psi?*?theta
????????????if?norm(x_r-x)<1e-6?%如果殘差小于1e-6則認(rèn)為恢復(fù)成功
????????????????P(8)?=?P(8)?+?1;
????????????end
???????end
???????for?iii?=?1:size(KIND,1)
???????????PercentageM(iii,kk)?=?P(iii)/CNT*100;?%計(jì)算恢復(fù)概率
???????end
????end
????for?jjj?=?1:size(KIND,1)
????????Percentage(1:length(K_set),mm,jjj)?=?PercentageM(jjj,:);
????end
end
toc
save?KtoPercentage1000gOMP?%運(yùn)行一次不容易，把變量全部存儲(chǔ)下來(lái)

%%?繪圖
S?=?['-ks';'-ko';'-yd';'-gv';'-b*';'-r.';'-rx';'-r+'];
figure;
for?mm?=?1:length(M_set)
????M?=?M_set(mm);
????K_set?=?5:5:70;
????L_Kset?=?length(K_set);
????for?ii?=?1:size(KIND,1)
????????plot(K_set,Percentage(1:L_Kset,mm,ii),S(ii,:));?%繪出x的恢復(fù)信號(hào)
????????hold?on;
????end
end
hold?off;
xlim([5?70]);
legend('OMP','ROMP','StOMP','SP','CoSaMP',...
????'gOMP(s=3)','gOMP(s=6)','gOMP(s=9)');
xlabel('Sparsity?level?K');
ylabel('The?Probability?of?Exact?Reconstruction');
title('Prob.?of?exact?recovery?vs.?the?signal?sparsity?K(M=128,N=256)(Gaussian)');

結(jié)論：gOMP只是在OMP基礎(chǔ)上修改了一下原子選擇的個(gè)數(shù)，效果就好很多。