




版權說明:本文檔由用戶提供并上傳,收益歸屬內容提供方,若內容存在侵權,請進行舉報或認領
文檔簡介
1、壓縮感知重構算法之正則化正交匹配追蹤(ROMP)正交匹配追蹤算法每次迭代均只選擇與殘差最相關的一列,自然人們會想:“每次迭代是否可以多選幾列呢”,正則化正交匹配追蹤(RegularizedOMP)就是其中一種改進方法。本篇將在上一篇壓縮感知重構算法之正交 匹配追蹤(OMP)的基礎上給出正則化正交匹配追蹤(ROMP)算法的MATLAB函數(shù)代碼,并且給出單次測試例程 代碼、測量數(shù)M與重構成功概率關系曲線繪制例程代碼。0、符號說明如下:壓縮觀測丫=血,其中j為觀測所得向量Mx1,*為原信號Nx1(M1e-6%判斷 productdes 中非零值個數(shù)break;endend%Identify:Choo
2、sea set J of the K biggest coordinatesif ii=KinJ = indexproductdes(1:Kin);%集合 JJval = productdes(1:Kin);%集合 J 對應的序列值K=Kin;else%or all of its nonzero coordinates,whichever is smaller TOC o 1-5 h z J=indexproductdes(1:ii);%集合JJval = productdes(1:ii);%集合J 對應的序列值K=ii;end%Regularize:AmongallsubsetsJ0 e J
3、 with comparable coordinatesMaxE =-1;%循環(huán)過程中存儲最大能量值for kk= 1:KJ0_tmp = zeros(1,K);iJ0 = 1;J0_tmp(iJ0) = J(kk);%以J(kk)為本次尋找J0的基準(最大值)Energy = Jval(kk)A2;%本次尋找 J0 的能量for mm = kk+1:Kif Jval(kk)2*Jval(mm)%找到符合|u(i)|=2|u(j)| 的iJ0 = iJ0 + 1;%J0 自變量增 1J0_tmp(iJ0) = J(mm);%更新 J0Energy = Energy + Jval(mm)A2;%
4、 更新能量else%不符合|u(i)|MaxE%本次所得J0的能量大于前一組J0 =J0_tmp(1:iJ0);%更新J0MaxE= Energy;%更新MaxE,為下次循環(huán)做準備endendpos=J0;val=productabs(J0);2、正則化正交匹配追蹤(ROMP)MATLAB代碼(CS_ROMP.m)這個函數(shù)是完全基于上一篇中的CS_OMP.m修改而成的。function theta = CS_ROMP( y,A,K )%CS_ROMP Summary of this function goes here%Version: 1.0 written by jbb0523 2015-
5、04-24% 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:Needell D, Vershynin R. Signal recovery from incomplete and% inaccurate measurements via regularized orthogonal matching pursuitJ.% IEEE Journal on Selected To
6、pics in Signal Processing, 2010, 4(2): 310316.y_rows,y_columns = size(y);if y_rows=2K)for ii=1:K%迭代 K 次product = A*r_n;%傳感矩陣A各列與殘差的內積%val,pos = max(abs(product);%找到最大內積絕對值,即與殘差最相關的列val,pos = Regularize(product,K);% 按正則化規(guī)則選擇原子At(:,Index+1:Index+length(pos) = A(:,pos);% 存儲這幾列Pos_theta(Index+1:Index+le
7、ngth(pos) = pos;% 存儲這幾列的序號if Index+length(pos)=M%At的行數(shù)大于列數(shù),此為最小二乘的基礎(列線性無關)Index = Index+length(pos);%更新 Index,為下次循環(huán)做準備else%At的列數(shù)大于行數(shù),列必為線性相關的,At(:,1:Index)*At(:,1:Index)將不可逆break;%跳出 for 循環(huán)endA(:,pos) = zeros(M,length(pos);%清零A的這幾列(其實此行可以不要因為它們與殘差正交)%y=At(:,1:Index)*theta,以下求 theta 的最小二乘解(Least Squ
8、are)theta_ls = (At(:,1:Index)*At(:,1:Index)A(-1)*At(:,1:Index)*y;% 最小二乘解%At(:,1:Index)*theta_ls 是 y 在 At(:,1:Index)列空間上的正交投影r_n = y - At(:,1:Index)*theta_ls;% 更新殘差if norm(r_n)=2*K%or until |I|=2K44.break;%44.break;%跳出for循環(huán)end45.endendtheta(Pos_theta(1:Index)=theta_ls;%恢復出的 thetaend3、ROMP單次重構測試代碼以下測試
9、代碼與上一篇中的OMP單次測試代碼基本完全一致,除了 M和K參數(shù)設置及調用CS_ROMP函數(shù)三處 之外。%壓縮感知重構算法測試clear all;close all;clc;M = 128;%觀測值個數(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*thetaPhi = randn(M,N);%測量矩陣為高斯矩陣A = Phi * Psi;
10、%傳感矩陣y = Phi * x;%得到觀測向量y%恢復重構信號xtictheta = CS_ROMP(y,A,K);x_r = Psi * theta;% x=Psi * thetatoc%繪圖figure;plot(x_r,k.-);%繪出x的恢復信號hold on;plot(x,r);%繪出原信號 xhold off;legend(Recovery,Original)fprintf(n 恢復殘差:);norm(x_r-x)%恢復殘差運行結果如下:(信號為隨機生成,所以每次結果均不一樣)1)圖:1510-550-101510-550-102) CommandwindowsElapsed t
11、ime (運行時間)is 0.022132 seconds.恢復殘差:ans=7.8066e-0154、測量數(shù)M與重構成功概率關系曲線繪制例程代碼以下測試代碼與上一篇中的OMP測量數(shù)M與重構成功概率關系曲線繪制例程代碼基本完全一致。本程序在 循環(huán)中填加了 “kk ”一行代碼并將“M = M_set(mm) ”一行的分號去掉,這是為了在運行過程中可以觀察程序運行狀 態(tài)、知道程序到哪一個位置。重構調用CS_ROMP函數(shù)并將save命令的文件名改為 ROMPMtoPercentage1000,以防止將OMP存儲的數(shù)據(jù)覆蓋(因為本人的所有代碼都在一個目錄下)。clear all;close all;c
12、lc;%參數(shù)配置初始化CNT = 1000;%對于每組(K,M,N),重復迭代次數(shù)N = 256;%信號x的長度Psi = eye(N);%x本身是稀疏的,定義稀疏矩陣為單位陣x=Psi*thetaK_set = 4,12,20,28,36;%信號 x 的稀疏度集合Percentage = zeros(length(K_set),N);% 存儲恢復成功概率%主循環(huán),遍歷每組(K,M,N)ticfor kk = 1:length(K_set)K = K_set(kk);% 本次稀疏度M_set = K:5:N;%M沒必要全部遍歷,每隔5測試一個就可以了PercentageK = zeros(1,
13、length(M_set);%存儲此稀疏度K下不同M的恢復成功概率kkfor mm = 1:length(M_set)M = M_set(mm)%本次觀測值個數(shù)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);%測量矩陣為高斯矩陣A=Phi * Psi;% 傳感矩陣y=Phi * x;%得到觀測向量ytheta = CS_ROMP(y,A,K);%恢復重構信號 thetax_
14、r = Psi * theta;% x=Psi * thetaif norm(x_r-x)1e-6%如果殘差小于1e-6則認為恢復成功P = P + 1;endendPercentageK(mm) = P/CNT*100;% 計算恢復概率endPercentage(kk,1:length(M_set) = PercentageK;endtocsave ROMPMtoPercentage1000 %運行一次不容易,把變量全部存儲下來%繪圖S = -ks;-ko;-kd;-kv;-k*;figure;for kk = 1:length(K_set)K = K_set(kk);M_set = K:5
15、:N;L_Mset = length(M_set);plot(M_set,Percentage(kk,1:L_Mset),S(kk,:);%繪出 x 的恢復信號hold on;endhold 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);本程序在聯(lián)想ThinkPadE430C筆記本(
16、4GBDDR3內存,i5-3210)上運行共耗時871.829395秒(上篇中OMP運行共耗時1171.606254秒),程序中將所有數(shù)據(jù)均通“save ROMPMtoPercentage1000存儲了下來, 以后可以再對數(shù)據(jù)進行分析,只需“l(fā)oad ROMPMtoPercentage1000”即可。程序運行結果比文獻4的Fig.1要差(OMP仿真結果比文獻中的要好),原因不詳。我調用了文獻2中的ROMP函數(shù),效果和這里的CS ROMP差不多。本程序運行結果;HEaxooaCDEElLICDECDdHEaxooaCDEElLICDECDd文獻4中的Fig.1:Percentage 寸 站胡si
17、gnals rqver-eciNurnber ofN)Percentage 寸 站胡signals rqver-eciNurnber ofN)Fig. 1 The percentage cf spar泥 fiat signals exactly rectjvtjred by ROMP as a fiinctiDn of the number of meaurmtnts N in dimensLon d = 256 for various levels of s-pusi ty np$#tl 必。Elc5、結語國內引用ROMP文獻時一般為4和5,文獻4更早一些,其實它們有關ROMP的敘述基本是一
18、樣子的, 下面分別是文獻4和5中的ROMP步驟:文獻4:Regularized Orthogonal Matching Pursuit (ROMP)Input: Measurement vector x E and sparsity level nOutput: Index set I C L ,dInitialized Let the index set 7 = 0 and the Tesidual r jc.Repeat the following steps until r = 0;Identify: Choose a set J of the n biggest coordijiate
19、s in magiutude of the observation vector u =如*廣,or all of its nonzero coordinates, whichever set is smaller.Regularize: Among all subsets J()C J with comparable coordinates2w(0| 2|m(j)| for all i, j 如choose 島 with the maximal energy l.Update: Add the set Jo to the index set: 1 JU 如 and update the re
20、sidual:y = arg min |x z|c; =工一y.zgRj文獻5:其實它們基本是一樣的,除了綜上迭代的條件略有不同,本文的CS_ROMP中都進行了考慮。其實若要看ROMP的原理,建議看文獻4的1.4節(jié):When we are trying to recover the signal u from its measurements x = we can use the observation vector u = *1: as a good local approximation to the signal v. Namely, the observation vector u e
21、ncodes correlations of the measurement vector x with the columns of .Note that 0 is a dictionary, and so since the signal v sparse, x has a sparse representation with respect to the dictionary. By the restricted isometry condition, every n columns form approximately an orthonormn system. Therefore,
22、every n coordinates of th巳 observation vector u look like correlations of the measurement vector x with the orthonormal basis audq therefore, are close in the Euclidean norm to the corresponding n coefficients of u. This is documented in Proposition 3.2 below.The local approximation property sugge&t
23、s to make us-e af the n biggest coordinates of the observation vector w, rather than one biggest coordinate as OMP did. We thus force the selected coordinates to be more regular (i.巳、closer to uniform) by selecting only the coordinates with comparable sizes. To this end, a new regufarizatim step will be needed lo ens.ure that each of these coordinates gets an even share of information. This leads to the following algorithm for sparse recovery:另外,針對ROMP的改進算法也很多,后面看些文獻再敘吧。推薦看看文獻67兩篇重構算法的綜述性文 獻,本文第一句話就來自文獻6。比較本文的ROMP仿真結果與上一篇中OMP的仿真結果可以發(fā)現(xiàn)效果遠沒有OMP好,當然本文文獻4 中的Fig.1比上一篇文獻1中的Fig.1,其中原因不詳。參考文獻
溫馨提示
- 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
- 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權益歸上傳用戶所有。
- 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內容里面會有圖紙預覽,若沒有圖紙預覽就沒有圖紙。
- 4. 未經(jīng)權益所有人同意不得將文件中的內容挪作商業(yè)或盈利用途。
- 5. 人人文庫網(wǎng)僅提供信息存儲空間,僅對用戶上傳內容的表現(xiàn)方式做保護處理,對用戶上傳分享的文檔內容本身不做任何修改或編輯,并不能對任何下載內容負責。
- 6. 下載文件中如有侵權或不適當內容,請與我們聯(lián)系,我們立即糾正。
- 7. 本站不保證下載資源的準確性、安全性和完整性, 同時也不承擔用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。
最新文檔
- 新疆師范大學《實驗室安全與法規(guī)》2023-2024學年第二學期期末試卷
- 重慶市主城區(qū)七校聯(lián)考2025年高三教學質量檢測試題試卷(二)物理試題含解析
- 公共交通運營服務收費標準制度
- 第18課 清朝的邊疆治理 教案2024-2025學年七年級歷史下冊新課標
- 內圓形吊頂施工方案
- 護坡植草施工方案
- 路基修復夜間施工方案
- 工程資料與施工方案
- 汽車隔音施工方案范本
- 2025年搞笑考試面試試題及答案
- 江蘇省藥品上市許可持有人藥品生產(chǎn)質量安全主體責任正面清單、負面清單(2023年版)
- 2024年GINA哮喘防治指南修訂解讀課件
- 木地板合同范本
- 2024中交二航局分包合同范本
- 2024年社區(qū)工作者考試必背1000題題庫必背(必刷)
- 教育改革與發(fā)展
- 《形體訓練》課件-勾繃腳訓練
- 醫(yī)療器械(耗材)項目投標服務實施投標方案(技術方案)
- 監(jiān)控系統(tǒng)維護保養(yǎng)方案
- 2023年國家廣播電視總局無線電臺管理局考試真題及答案
- 房屋修繕工程技術規(guī)程 DG-TJ08-207-2008
評論
0/150
提交評論