基于壓縮感知的正交匹配算法圖像重建摘要：壓縮感知理論是由Donoho和Candes提出的一種充分利用信號稀疏性的全新的信號采樣理論。該理論表明，用遠低于Nyquist采樣定理要求的頻率對信號進行采樣也能實現信號的精確重構。該理論突破了傳統的以Nyquist定理為基準的信號處理方法，實現了在獲取數據的同時對其進行適當的壓縮，克服了采樣數據量大，采樣時間長及數據存儲空間浪費嚴重的問題，因此進一步降低了信號處理的時間和器件成本。壓縮感知理論有三個核心方面：（1）稀疏變換，即對一個非稀疏的信號，找到一個合適的正交基使該信號在它上可以稀疏表示；（2）測量矩陣，與變換基不相干且平穩的矩陣；（3）重構算法，利用數學算法完成對信號的精確重構，該過程可看為求解一個優化問題。本文介紹了主要介紹了壓縮感知原理和目前最為成熟的壓縮感知重建算法——正交匹配追蹤算法，通過MATLAB平臺設計實現了基本的正交匹配追蹤算法，對一維、二維信號進行了重建仿真。關鍵詞：壓縮感知；稀疏變換；正交匹配；圖像重建BasedOnCompressedSensingOfOrthogonalMatchingAlgorithmImageRecoveryAbstract：CompressedsensingisanovelsamplingtheorywhichisproposedbyDonohoandCandès.Thistheoryisundertheconditionthatthesignaliscompressibleorsparse.Inthiscase,usingfarlessthantherequiredsamplingfrequencyoftheNyquisttheorytosamplethesignalisabletoaccuratelyreconstructthesignal.CompressedtheorybreaksthoughthetraditionalNyquistsamplingtheory,whichovercomesalotofproblemssuchasagreatnumberofsamplingdata,timewasting,datastoragespacewastingandsoon.Asaresult,itreducessignalprocessingcostanddevicecost.Thecompressedtheoryhasthreekeysides:(1)Sparsetransformation,foranon-sparsesignal,weneedtofindaproperorthogonalbasisonwhichthesignalhasasparserepresentation;(2)Observationmatrix,itisirrelevantwiththeorthogonalbasis;(3)reconstructionalgorithms,usingareconstructionalgorithmtoensuretheaccuracyofthesignalreconstruction,thewholeprocesscanbeconsideredasthesolvetoaoptimizationproblem.ThispaperintroducesCSandmostmaturecompressionperceptionalgorithmatpresent-Orthogonalmatchingalgorithm.ThroughtheMATLABdesignrealizebasicorthogonalmatchingalgorithms,ThroughtheMATLABdesignrealizebasicorthogonalmatchingalgorithmofone-dimensional,two-dimensionalsignalprocessingsimulation.Keywords：Compressedsensing;Sparsetransform;Orthogonalmatching;Imagerecovery.西安文理學院本科畢業設計（論文）第頁參考文獻[1]李樹濤,魏丹.壓縮傳感綜述.自動化學報,2009,35(11):1369-1377.[2]邵文澤,韋志輝,肖亮等.壓縮感知基本理論：回顧與展望.中國科技論文在線.[3]金堅,谷源濤,梅順良.壓縮采樣技術及應用.電子與信息學報,2010,62(2):470-475[4]石光明.劉丹華.高大化.劉哲.林杰.王良君壓縮感知理論及其研究進展-ACTAElectronicaSinica2009，37(5)[5]沙威.“壓縮傳感”引論.http://www.eee.hku.hk/~wsha[6]何雪云,宋榮芳,周克琴.基于壓縮感知的OFDM系統系數信道估計新方法研究.南京郵電大學學報,2010,30(2):60-65[7]JustinRomberg.ImagingviaCompressiveSampling.IEEESignalProcessingMagazine,2008,14-20[8]張銳基于壓縮感知理論的圖像壓縮初步研究-ComputerKnowledgeAndTechnology2010，6(4)[9]喻玲娟.謝曉春壓縮感知理論簡介-VideoEngineering2008，32(12)[10]張春梅.尹忠科.肖明霞基于冗余字典的信號超完備表示與稀疏分解-科學通報2006(06)[11]JoelA.TroppandAnnaC.GilbertSignalRecoveryFromRandomMeasurementsViaOrthogonalMatchingPursuit，IEEETRANSACTIONSONINFORMATIONTHEORY,VOL.53,NO.12[12]DONOHOD.TSAIGYExtensionsofcompressedsensing2006(03)[13]ECandesandJRomberg,Quantitativerobustuncentaintyprinciplesandoptimallysparsedecompositions[J].FoundationsofComputMath,2006,6(2):227-254．[14]李小波.基于壓縮感知的測量矩陣研究.北京交通大學碩士學位論文.2010:6[15]石光明，劉丹華，高大化.壓縮感知理論及研究進展.電子學報.2009,37(5):1075-1076附錄一源程序清單程序1：一維原始信號的生成程序：clc;clear%%1.時域測試信號生成K=7;%稀疏度(做FFT可以看出來)N=256;%信號長度M=64;%測量數(M>=K*log(N/K),至少40,但有出錯的概率)f1=50;%信號頻率1f2=100;%信號頻率2f3=200;%信號頻率3f4=400;%信號頻率4fs=800;%采樣頻率ts=1/fs;%采樣間隔Ts=1:N;%采樣序列x=0.3*sin(2*pi*f1*Ts*ts)+0.6*sin(2*pi*f2*Ts*ts)+0.1*sin(2*pi*f3*Ts*ts)+0.9*sin(2*pi*f4*Ts*ts);%完整信號plot(x,'r')%原始信號程序2：一維重建信號的生成clc;clear%%1.時域測試信號生成K=7;%稀疏度(做FFT可以看出來)N=256;%信號長度M=64;%測量數(M>=K*log(N/K),至少40,但有出錯的概率)f1=50;%信號頻率1f2=100;%信號頻率2f3=200;%信號頻率3f4=400;%信號頻率4fs=800;%采樣頻率ts=1/fs;%采樣間隔Ts=1:N;%采樣序列x=0.3*sin(2*pi*f1*Ts*ts)+0.6*sin(2*pi*f2*Ts*ts)+0.1*sin(2*pi*f3*Ts*ts)+0.9*sin(2*pi*f4*Ts*ts);%完整信號%%2.時域信號壓縮傳感Phi=randn(M,N);%測量矩陣(高斯分布白噪聲)s=Phi*x.';%獲得線性測量%%3.正交匹配追蹤法重構信號(本質上是1-范數最優化問題)m=2*K;%算法迭代次數(m>=K)Psi=fft(eye(N,N))/sqrt(N);%傅里葉正變換矩陣T=Phi*Psi';%恢復矩陣(測量矩陣*正交反變換矩陣)hat_y=zeros(1,N);%待重構的譜域(變換域)向量Aug_t=[];%增量矩陣(初始值為空矩陣)r_n=s;%殘差值fortimes=1:m;%迭代次數(有噪聲的情況下,該迭代次數為K)forcol=1:N;%恢復矩陣的所有列向量product(col)=abs(T(:,col)'*r_n);%恢復矩陣的列向量和殘差的投影系數(內積值)end[val,pos]=max(product);%最大投影系數對應的位置Aug_t=[Aug_t,T(:,pos)];%矩陣擴充T(:,pos)=zeros(M,1);%選中的列置零（實質上應該去掉，為了簡單我把它置零）aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;%最小二乘,使殘差最小r_n=s-Aug_t*aug_y;%殘差pos_array(times)=pos;%紀錄最大投影系數的位置endhat_y(pos_array)=aug_y;%重構的譜域向量hat_x=real(Psi'*hat_y.');%做逆傅里葉變換重構得到時域信號%%4.恢復信號和原始信號對比figure(1)holdon;plot(hat_x,'k.-')%重建信號plot(x,'r')%原始信號legend('Recovery','Original')norm(hat_x.'-x)/norm(x)%重構誤差程序3：二維圖像OMP重建functionWavelet_OMPclc;clearX=imread('lena256.bmp');%讀文件X=double(X);[a,b]=size(X);ww=DWT(a);%小波變換矩陣生成X1=ww*sparse(X)*ww';%小波變換讓圖像稀疏化（注意該步驟會耗費時間，但是會增大稀疏度X1=full(X1);M=190;%隨機矩陣生成R=randn(M,a);Y=R*X1;%測量%OMP算法X2=zeros(a,b);%恢復矩陣fori=1:b%列循環rec=omp(Y(:,i),R,a);X2(:,i)=rec;endfigure(1);%原始圖像imshow(uint8(X));title('原始圖像');figure(2);%變換圖像imshow(uint8(X1));title('小波變換后的圖像');figure(3);%壓縮傳感恢復的圖像X3=ww'*sparse(X2)*ww;%小波反變換X3=full(X3);imshow(uint8(X3));title('恢復的圖像');%誤差(PSNR)errorx=sum(sum(abs(X3-X).^2));%MSE誤差psnr=10*log10(255*255/(errorx/a/b))%PSNR%OMP的函數%s-測量；T-觀測矩陣；N-向量大小functionhat_y=omp(s,T,N)Size=size(T);%觀測矩陣大小M=Size(1);%測量hat_y=zeros(1,N);%待重構的譜域(變換域)向量Aug_t=[];%增量矩陣(初始值為空矩陣)r_n=s;%殘差值fortimes=1:M/4;%迭代次數(稀疏度是測量的1/4)forcol=1:N;%恢復矩陣的所有列向量product(col)=abs(T(:,col)'*r_n);%恢復矩陣的列向量和殘差的投影系數(內積值)end[val,pos]=max(product);%最大投影系數對應的位置Aug_t=[Aug_t,T(:,pos)];%矩陣擴充T(:,pos)=zeros(M,1);%選中的列置零（實質上應該去掉，為了簡單我把它置零）aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;%最小二乘,使殘差最小r_n=s-Aug_t*aug_y;%殘差pos_array(times)=pos;%紀錄最大投影系數的位置if(norm(r_n)<9)%殘差足夠小break;endendhat_y(pos_array)=aug_y;%重構的向量程序4：BP、OMP、STOMP_FDR重建圖像x=imread('lena.bmp');%讀文件[m,n]=size(x);xrec_BP=size(x);xrec_OMP=size(x);xrec_FDR=size(x);t_BP=0;t_OMP=0;t_FDR=0;fori=1:mx1=x(i,:);[cA,cD]=dwt(x1,'db1');%小波變換c1=length(cD)Mdetail=c1;Ndetail=floor(c1*0.4);A=randn(Ndetail,Mdetail);%隨機矩陣生成y=A*cD';q=0.5;S=10;tic;alpha=SolveBP(A,y,Mdetail);%BP算法處理圖像t_BP=t_BP+toc;rec1=idwt(cA,alpha','db1');xrec_BP(i,1:n)=rec1;tic;[alpha,iters,activeset]=SolveOMP(A,y,Mdetail);%OMP算法處理圖像t_OMP=t_OMP+toc;rec2=idwt(cA,alpha','db1');xrec_OMP(i,1:n)=rec2;tic;[alpha,iters]=SolveStOMP(A,y,Mdetail,'FDR',q,S);%STOMP_FDR算法處理圖像t_FDR=t_FDR+toc;rec3=idwt(cA,alpha','db1');xrec_FDR(i,1:n)=rec3;endsubplot(2,2,1);imshow(x);title(['Origineimage','N=',num2str(m*n)]);%原始圖像subplot(2,2,2);imshow(uint8(xrec_BP));title(['BP,','samp=',num2str((Ndetail/Mdetail)*100),'%time=',num2str(t_BP),'sec']);%BP恢復圖像subplot(2,2,3);imshow(uint8(xrec_OMP));title(['OMP,','samp=',num2str((Ndetail/Mdetail)*100),'%time=',num2str(t_OMP),'sec']);%OMP恢復圖像subplot(2,2,4);imshow(uint8(xrec_FDR));title(['FDR,','samp=',num2str((Ndetail/Mdetail)*100),'%time=',num2str(t_FDR),'sec']);%STOMP恢復圖像附錄二英文文獻翻譯英文文獻：ImageProcessingImageprocessingisnotaonestepprocess．Weareabletodistinguishbetweenseveralstepswhichmustbeperformedoneaftertheotheruntilwecanextractthedataofinterestfromtheobservedscene．Imageprocessingbeginswiththecaptureofanimagewithasuitable，notnecessarilyoptical，acquisitionsystem．Inatechnicalorscientificapplication，wemaychoosetoselectanappropriateimagingsystem．Furthermore，wecansetuptheilluminationsystem，choosethebestwavelengthrange，andselectotheroptionstocapturetheobjectfeatureofinterestinthebestwayinanimage．Oncetheimageissensed，itmustbebroughtintoaformthatcanbetreatedwithdigitalcomputers.Thisprocessiscalleddigitization.Thefirststepsofdigitalprocessingmayincludeanumberofdifferentoperationsandareknownasimageprocessing．Ifthesensorhasnonlinearcharacteristics,these

needtobecorrected．Likewise，brightnessandcontrastoftheimagemayrequireimprovement．Commonly，too，coordinatetransformationsareneededtorestoregeometricaldistortionsintroducedduringimageformation．Radiometricandgeometriccorrectionsareelementarypixelprocessingoperations．Itmaybenecessarytocorrectknowndisturbancesintheimage，forinstancecausedbyadefocusedoptics，motionblur，errorsinthesensor，orerrorsinthetransmissionofimagesignals．Wealsodealwithreconstructiontechniqueswhicharerequiredwithmanyindirectimagingtechniquessuchastomographythatdelivernodirectimage.Awholechainofprocessingstepsisnecessarytoanalyzeandidentifyobjects．First，adequatefilteringproceduresmustbeappliedinordertodistinguishtheobjectsofinterestfromotherobjectsandthebackground．Essentially，fromanimage（orseveralimages），oneormorefeatureimagesareextracted．Thebasictoolsforthistaskareaveragingandedgedetectionandtheanalysisofsimpleneighborhoodsandcomplexpatternsknownastextureinimageprocessing．Animportantfeatureofanobjectisalsoitsmotion．Techniquestodetectanddeterminemotionarenecessary．Thentheobjecthastobeseparatedfromthebackground．Thismeansthatregionsofconstantfeaturesanddiscontinuitiesmustbeidentified．Thisprocessleadstoalabelimage．Nowthatweknowtheexactgeometricalshapeoftheobject，wecanextractfurtherinformationsuchasthemeangrayvalue，thearea，perimeter，andotherparametersfortheformoftheobject．Theseparameterscanbeusedtoclassifyobjects．Thisisanimportantstepinmanyapplicationsofimageprocessing，asthefollowingexamplesshow：Inasatelliteimageshowinganagriculturalarea，wewouldliketodistinguishfieldswithdifferentfruitsandobtainparameterstoestimatetheirripenessortodetectdamagebyparasites．Therearemanymedicalapplicationswheretheessentialproblemistodetectpathologi-alchanges．Aclassicexampleistheanalysisofaberrationsinchromosomes．Characterrecognitioninprintedandhandwrittentextisanotherexamplewhichhasbeenstudiedsinceimageprocessingbeganandstillposessignificantdifficulties．Youhopefullydomore，namelytrytounderstandthemeaningofwhatyouarereading．Thisisalsothefinalstepofimageprocessing，whereoneaimstounderstandtheobservedscene．Weperformthistaskmoreorlessunconsciouslywheneverweuseourvisualsystem．Werecognizepeople，wecaneasilydistinguishbetweentheimageofascientificlabandthatofalivingroom，andwewatchthetraffictocrossastreetsafely．Wealldothiswithoutknowinghowthevisualsystemworks．Forsometimesnow，imageprocessingandcomputer-graphicshavebeentreatedastwodifferentareas．Knowledgeinbothareashasincreasedconsiderablyandmorecomplexproblemscannowbetreated．Computergraphicsisstrivingtoachievephotorealisticcomputer-generatedimagesofthree-dimensionalscenes，whileimageprocessingistryingtoreconstructonefromanimageactuallytakenwithacamera．Inthissense，imageprocessingperformstheinverseproceduretothatofcomputergraphics．Westartwithknowledgeoftheshapeandfeaturesofanobject.andworkupwardsuntilwegetatwo-dimensionalimage．Tohandleimageprocessingorcomputergraphics，webasicallyhavetoworkfromthesameknowledge．Weneedtoknowtheinteractionbetweenilluminationandobjects，howathree-dimensionalsceneisprojectedontoanimageplane，etc．Therearestillquiteafewdifferencesbetweenanimageprocessingandagraphicsworkstation．Butwecanenvisagethat，whenthesimilaritiesandinterrelationsbetweencomputergraphicsandimageprocessingarebetterunderstoodandtheproperhardwareisdeveloped，wewillseesomekindofgeneral-purposeworkstationinthefuturewhichcanhandlecomputergraphicsaswellasimageprocessingtasks．Theadventofmultimedia，I.e.，theintegrationoftext，images，sound，andmovies，willfurtheracceleratetheunificationofcomputergraphicsandimageprocessing.InJanuary1980ScientificAmericanpublishedaremarkableimagecalledPlume2，thesecondofeightvolcaniceruptionsdetectedontheJovianmoonbythespacecraftVoyager1on5March1979．Thepicturewasalandmarkimageininterplanetaryexploration—thefirsttimeaneruptingvolcanohadbeenseeninspace．Itwasalsoatriumphforimageprocessing．Satelliteimageryandimagesfrominterplanetaryexplorershaveuntilfairlyrecentlybeenthemajorusersofimageprocessingtechniques，whereacomputerimageisnumericallymanipulatedtoproducesomedesiredeffect-suchasmakingaparticularaspectorfeatureintheimagemorevisible.ImageprocessinghasitsrootsinphotoreconnaissanceintheSecondWorldWarwhereprocessingoperationswereopticalandinterpretationoperationswereperformedbyhumanswhoundertooksuchtasksasquantifyingtheeffectofbombingraids．Withtheadventofsatelliteimageryinthelate1960s，muchcomputer-basedworkbeganandthecolorcompositesatelliteimages，sometimesstartlinglybeautiful,havebecomepartofourvisualcultureandtheperceptionofourplanet．Likecomputergraphics，itwasuntilrecentlyconfinedtoresearchlaboratorieswhichcouldaffordtheexpensiveimageprocessingcomputersthatcouldcopewiththesubstantialprocessingoverheadsrequiredtoprocesslargenumbersofhigh-resolutionimages．Withtheadventofcheappowerfulcomputersandimagecollectiondeviceslikedigitalcamerasandscanners，wehaveseenamigrationofimageprocessingtechniquesintothepublicdomain．Classicalimageprocessingtechniquesareroutinelyemployedbygraphicdesignerstomanipulatephotographicandgeneratedimagery，eithertocorrectdefects，changecolorandsoonorcreativelytotransformtheentirelookofanimagebysubjectingittosomeoperationsuchasedgeenhancement.Arecentmainstreamapplicationofimageprocessingisthecompressionofimages—eitherfortransmissionacrosstheInternetorthecompressionofmovingvideoimagesinvideotelephonyandvideoconferencing．Videotelephonyisoneofthecurrentcrossoverareasthatemploybothcomputergraphicsandclassicalimageprocessingtechniquestotrytoachieveveryhighcompressionrates．Allthisispartofaninexorabletrendtowardsthedigitalrepresentationofimages．Indeedthatmostpowerfulimageformofthetwentiethcentury—theTVimage—isalsoabouttobetakenintothedigitaldomain．Imageprocessingischaracterizedbyalargenumberofalgorithmsthatarespecificsolutionstospecificproblems．Somearemathematicalorcontext-independentoperationsthatareappliedtoeachandeverypixel．Forexample，wecanuseFouriertransformstoperformimagefilteringoperations．Othersare“algorithmic”—wemayuseacomplicatedrecursivestrategytofindthosepixelsthatconstitutetheedgesinanimage．Imageprocessingoperationsoftenformpartofacomputervisionsystem．Theinputimagemaybefilteredtohighlightorrevealedgespriortoashapedetectionusuallyknownaslow-leveloperations．Incomputergraphicsfilteringoperationsareusedextensivelytoavoidabasingorsamplingartifacts.中文翻譯：圖像處理圖像處理不是一步就能完成的過程?？蓪⑺殖芍T多步驟，必須一個接一個地執行這些步驟，直到從被觀察的景物中提取出有用的數據。圖像處理首先是以適當的但不一定是光學的采集系統對圖像進行采集。在技術或科學應用中，可以選擇一個適當的成像系統。此外，可以建立照明系統，選擇最佳波長范圍，以及選擇其他方案以便用最好的方法在圖像中獲取有用的對象特征。一旦圖像被檢測到，必須將其變成數字計算機可處理的形式，這個過程稱之為數字化。數字化處理的第一步包含了一系列不同的操作并被稱之為圖像處理。如果傳感器具有非線性特性，就必須予以校正，同樣，圖像的亮度和對比度也需要改善。通常，還需要進行坐標變換以消除在成像時產生的幾何畸變。輻射度校正和幾何校正是最基本的像素處理操作。在圖像中，對已知的干擾進行校正也是不可少的，比如由于光學聚焦不準，運動模糊，傳感器誤差以及圖像信號傳輸誤差所引起的干擾。在此還要涉及圖像重構技術，它需要許多間接的成像技術，比如不直接提供圖像的X射線斷層技術等。一套完整的處理步驟對于物體的分析和識別是必不可少的。首先，應該采用適當的過濾技術以便從其他物體和背景中將所感興趣的物體區分出來。實質上就是從一幅圖像（或者數幅圖像）中抽取出一幅或幾幅特征圖像。要完成這個任務最基本的工具就是圖像處理中所使用的求均值和邊緣檢測、簡單的相鄰像素分析，以及復雜的被稱為材質描述的模式分析。物體的一個重要特性就是它的運動性。檢測和確定物體運動性的技術是必不可少的。隨后，該物體必須從背景中分離出來，這就意味著具有同樣特性和不同特性的區域必須被識別出來。這個過程產生出標志圖像。既然已經知道了物體精確的幾何形狀，就可以抽取諸如平均灰度值、區域、邊界以及形成物體的其他參數等更多的信息。這些參數可用來對物體進行分類，這是許多圖像處理應用中至關重要的一步，比如下面一些應用：在一個顯示農業