淺談壓縮感知（二十五）：壓縮感知重構(gòu)算法之分段正交匹配追蹤（StOMP）

主要內(nèi)容：

StOMP的算法流程StOMP的MATLAB實現(xiàn)一維信號的實驗與結(jié)果門限參數(shù)Ts、測量數(shù)M與重構(gòu)成功概率關(guān)系的實驗與結(jié)果一、StOMP的算法流程

分段正交匹配追蹤(Stagewise OMP)也是由OMP改進而來的一種貪心算法，與CoSaMP、SP算法類似，不同之處在于CoSaMP、SP算法在迭代過程中選擇的是與信號內(nèi)積最大的2K或K個原子，而StOMP是通過門限閾值來確定原子。此算法的輸入?yún)?shù)中沒有信號稀疏度K，因此相比于ROMP及CoSaMP有獨到的優(yōu)勢（這句話存在疑問）。

StOMP的算法流程：

二、StOMP的MATLAB實現(xiàn)(CS_StOMP.m)

function?[?theta?]?=?CS_StOMP(?y,A,S,ts?)
%???CS_StOMP
%???Detailed?explanation?goes?here
%???y?=?Phi?*?x
%???x?=?Psi?*?theta
%????y?=?Phi*Psi?*?theta
%???令?A?=?Phi*Psi,?則y=A*theta
%???S?is?the?maximum?number?of?StOMP?iterations?to?perform
%???ts?is?the?threshold?parameter
%???現(xiàn)在已知y和A，求theta
%???Reference:Donoho?D?L，Tsaig?Y，Drori?I，Starck?J?L．Sparse?solution?of
%???underdetermined?linear?equations?by?stagewise?orthogonal?matching?
%???pursuit[J]．IEEE?Transactions?on?Information?Theory，2012，58(2)：1094—1121
????if?nargin?<?4
????????ts?=?2.5;?%ts范圍[2,3],默認值為2.5
????end
????if?nargin?<?3
????????S?=?10;?%S默認值為10
????end
????[y_rows,y_columns]?=?size(y);
????if?y_rowsts*sigma);?%選出大于閾值的列
????????Is?=?union(pos_num,Js);?%pos_num與Js并集
????????if?length(pos_num)?==?length(Is)
????????????if?ss==1
????????????????theta_ls?=?0;?%防止第1次就跳出導致theta_ls無定義
????????????end
????????????break;?%如果沒有新的列被選中則跳出循環(huán)
????????end
????????%At的行數(shù)要大于列數(shù)，此為最小二乘的基礎(chǔ)(列線性無關(guān))
????????if?length(Is)<=M
????????????pos_num?=?Is;?%更新列序號集合
????????????At?=?A(:,pos_num);?%將A的這幾列組成矩陣At
????????else?%At的列數(shù)大于行數(shù)，列必為線性相關(guān)的,At'*At將不可逆
????????????if?ss==1
????????????????theta_ls?=?0;?%防止第1次就跳出導致theta_ls無定義
????????????end
????????????break;?%跳出for循環(huán)
????????end
????????%y=At*theta，以下求theta的最小二乘解(Least?Square)
????????theta_ls?=?(At'*At)^(-1)*At'*y;?%最小二乘解
????????%At*theta_ls是y在At列空間上的正交投影
????????res?=?y?-?At*theta_ls;?%更新殘差
????????if?norm(res)<1e-6?%Repeat?the?steps?until?r=0
????????????break;?%跳出for循環(huán)
????????end
????end
????theta(pos_num)=theta_ls;?%恢復(fù)出的theta
end

三、一維信號的實驗與結(jié)果

%壓縮感知重構(gòu)算法測試
clear?all;close?all;clc;
M?=?64;?%觀測值個數(shù)
N?=?256;?%信號x的長度
K?=?12;?%信號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)/sqrt(M);?%測量矩陣為高斯矩陣
A?=?Phi?*?Psi;?%傳感矩陣
y?=?Phi?*?x;?%得到觀測向量y

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

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

四、門限參數(shù)ts、測量數(shù)M與重構(gòu)成功概率關(guān)系的實驗與結(jié)果

clear?all;close?all;clc;

%%?參數(shù)配置初始化
CNT?=?1000;%對于每組(K,M,N)，重復(fù)迭代次數(shù)
N?=?256;%信號x的長度
Psi?=?eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta
ts_set?=?2:0.2:3;
K_set?=?[4,12,20,28,36];%信號x的稀疏度集合
Percentage?=?zeros(N,length(K_set),length(ts_set));%存儲恢復(fù)成功概率

%%?主循環(huán)，遍歷每組(ts,K,M,N)
tic
for?tt?=?1:length(ts_set)
????ts?=?ts_set(tt);
????for?kk?=?1:length(K_set)
????????K?=?K_set(kk);%本次稀疏度
????????%M沒必要全部遍歷，每隔5測試一個就可以了
????????M_set=2*K:5:N;
????????PercentageK?=?zeros(1,length(M_set));%存儲此稀疏度K下不同M的恢復(fù)成功概率
????????for?mm?=?1:length(M_set)
???????????M?=?M_set(mm);%本次觀測值個數(shù)
???????????fprintf('ts=%f,K=%d,M=%dn',ts,K,M);
???????????P?=?0;
???????????for?cnt?=?1:CNT?%每個觀測值個數(shù)均運行CNT次
????????????????Index_K?=?randperm(N);
????????????????x?=?zeros(N,1);
????????????????x(Index_K(1:K))?=?5*randn(K,1);%x為K稀疏的，且位置是隨機的????????????????
????????????????Phi?=?randn(M,N)/sqrt(M);%測量矩陣為高斯矩陣
????????????????A?=?Phi?*?Psi;%傳感矩陣
????????????????y?=?Phi?*?x;%得到觀測向量y
????????????????theta?=?CS_StOMP(y,A,10,ts);%恢復(fù)重構(gòu)信號theta
????????????????x_r?=?Psi?*?theta;%?x=Psi?*?theta
????????????????if?norm(x_r-x)<1e-6%如果殘差小于1e-6則認為恢復(fù)成功
????????????????????P?=?P?+?1;
????????????????end
???????????end
???????????PercentageK(mm)?=?P/CNT*100;%計算恢復(fù)概率
????????end
????????Percentage(1:length(M_set),kk,tt)?=?PercentageK;
????end
end
toc
save?StOMPMtoPercentage1000?%運行一次不容易，把變量全部存儲下來

%%?繪圖
for?tt?=?1:length(ts_set)
????S?=?['-ks';'-ko';'-kd';'-kv';'-k*'];
????figure;
????for?kk?=?1:length(K_set)
????????K?=?K_set(kk);
????????M_set=2*K:5:N;
????????L_Mset?=?length(M_set);
????????plot(M_set,Percentage(1:L_Mset,kk,tt),S(kk,:));%繪出x的恢復(fù)信號
????????hold?on;
????end
????hold?off;
????xlim([0?256]);
????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,ts=',...
????????num2str(ts_set(tt)),')(Gaussian)']);
end
for?kk?=?1:length(K_set)
????K?=?K_set(kk);
????M_set=2*K:5:N;
????L_Mset?=?length(M_set);
????S?=?['-ks';'-ko';'-kd';'-kv';'-k*';'-k+'];
????figure;
????for?tt?=?1:length(ts_set)
????????plot(M_set,Percentage(1:L_Mset,kk,tt),S(tt,:));%繪出x的恢復(fù)信號
????????hold?on;
????end
????hold?off;
????xlim([0?256]);
????legend('ts=2.0','ts=2.2','ts=2.4','ts=2.6','ts=2.8','ts=3.0');
????xlabel('Number?of?measurements(M)');
????ylabel('Percentage?recovered');
????title(['Percentage?of?input?signals?recovered?correctly(N=256,K=',...
????????num2str(K),')(Gaussian)']);????
end

1、門限參數(shù)ts分別為2.0,2.2,2.4,2.6,2.8,3.0時，不同稀疏信號下，測量值M與重構(gòu)成功概率的關(guān)系：

2、稀疏度為4,12,20,28,36時，不同門限參數(shù)ts下，測量值M與重構(gòu)成功概率的關(guān)系：

??

結(jié)論：

通過對比可以看出，總體上講ts=2.4或ts=2.6時效果較好，較大和較小重構(gòu)效果都會降低，這里由于沒有ts=2.5的情況，但我們推測ts=2.5應(yīng)該是一個比較好的值，因此一般默認取為2.5即可。