主要內容：
- CoSaMP的算法流程
- CoSaMP的MATLAB實現
- 一維信號的實驗與結果
- 測量數M與重構成功概率關系的實驗與結果
一、CoSaMP的算法流程
壓縮采樣匹配追蹤(CompressiveSampling MP)是D. Needell繼ROMP之后提出的又一個具有較大影響力的重構算法。CoSaMP也是對OMP的一種改進，每次迭代選擇多個原子，除了原子的選擇標准之外，它有一點不同於ROMP：ROMP每次迭代已經選擇的原子會一直保留，而CoSaMP每次迭代選擇的原子在下次迭代中可能會被拋棄。
二、CS_CoSaMP的MATLAB實現(CS_CoSaMP.m)
function [ theta ] = CS_CoSaMP( y,A,K ) % CS_CoSaOMP % Detailed explanation goes here % y = Phi * x % x = Psi * theta % y = Phi*Psi * theta % 令 A = Phi*Psi, 則y=A*theta % K is the sparsity level % 現在已知y和A，求theta % Reference:Needell D，Tropp J A．CoSaMP：Iterative signal recovery from % incomplete and inaccurate samples[J]．Applied and Computation Harmonic % Analysis，2009，26：301-321. [m,n] = size(y); if m<n y = y'; %y should be a column vector end [M,N] = size(A); %傳感矩陣A為M*N矩陣 theta = zeros(N,1); %用來存儲恢復的theta(列向量) pos_num = []; %用來迭代過程中存儲A被選擇的列序號 res = y; %初始化殘差(residual)為y for kk=1:K %最多迭代K次 %(1) Identification product = A'*res; %傳感矩陣A各列與殘差的內積 [val,pos]=sort(abs(product),'descend'); Js = pos(1:2*K); %選出內積值最大的2K列 %(2) Support Merger Is = union(pos_num,Js); %Pos_theta與Js並集 %(3) Estimation %At的行數要大於列數，此為最小二乘的基礎(列線性無關) if length(Is)<=M At = A(:,Is); %將A的這幾列組成矩陣At else %At的列數大於行數，列必為線性相關的,At'*At將不可逆 if kk == 1 theta_ls = 0; end break; %跳出for循環 end %y=At*theta，以下求theta的最小二乘解(Least Square) theta_ls = (At'*At)^(-1)*At'*y; %最小二乘解 %(4) Pruning [val,pos]=sort(abs(theta_ls),'descend'); %(5) Sample Update pos_num = Is(pos(1:K)); theta_ls = theta_ls(pos(1:K)); %At(:,pos(1:K))*theta_ls是y在At(:,pos(1:K))列空間上的正交投影 res = y - At(:,pos(1:K))*theta_ls; %更新殘差 if norm(res)<1e-6 %Repeat the steps until r=0 break; %跳出for循環 end end theta(pos_num)=theta_ls; %恢復出的theta end
三、一維信號的實驗與結果
%壓縮感知重構算法測試 clear all;close all;clc; M = 64; %觀測值個數 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); %測量矩陣為高斯矩陣 A = Phi * Psi; %傳感矩陣 y = Phi * x; %得到觀測向量y %% 恢復重構信號x tic theta = CS_CoSaMP( 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) %恢復殘差
四、測量數M與重構成功概率關系的實驗與結果
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 for kk = 1:length(K_set) K = K_set(kk); %本次稀疏度 M_set = 2*K:5:N; %M沒必要全部遍歷，每隔5測試一個就可以了 PercentageK = zeros(1,length(M_set)); %存儲此稀疏度K下不同M的恢復成功概率 for mm = 1:length(M_set) M = M_set(mm); %本次觀測值個數 fprintf('K=%d,M=%d\n',K,M); P = 0; for cnt = 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)/sqrt(M); %測量矩陣為高斯矩陣 A = Phi * Psi; %傳感矩陣 y = Phi * x; %得到觀測向量y theta = CS_CoSaMP(y,A,K); %恢復重構信號theta x_r = Psi * theta; % x=Psi * theta if norm(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 CoSaMPMtoPercentage1000 %運行一次不容易，把變量全部存儲下來 %% 繪圖 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(kk,1:L_Mset),S(kk,:));%繪出x的恢復信號 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)(Gaussian)');
