压缩感知入门matlab代码

BeginnersCodeforCompressiveSensingAlejandroWeinsteinSeptember20091SparseSignalsintheTimeDomain1.1UsingaRandomSensingMatrixInthis¯rstexamplewewillmeasureasignalthatissparseinthetimedomain.Wewillusearandomsensingmatrix,andwewillsolvetherecoveryproblemusingthel1-Magictoolbox.Weusethefollowingfunctionstogeneratethesignalsandthesensingmatrix:Listing1:Sparsesignalandrandommeasurementmatrix.functionf=getsparsefun(n,s)tmp=randperm(n);f=zeros(n,1);f(tmp(1:s))=randn(s,1);functionA=getArandom(n,m)A=sqrt(1/m)¤randn(m,n);Thefollowingscriptusethesefunctionstogeneratethesignal,takethemeasurementsanddotherecovery.Figure1showstheresult.Listing2:Example1.1%CSexample12%Sensingmatrixphiisrandom.3%RepresentationbasisPsiisthecanonicalbasis.4%Recoveringusingl1magic.56clc7clearall8closeall910n=512;%Signallength11s=25;%Sparsitylevel12m=5¤s;%Numberofmeasurements1314f=getsparsefun(n,s);15A=getArandom(n,m);1617y=A¤f;%Takethemeasurements1819%Solveusingl1magic.20path(path,'./Optimization');21x0=pinv(A)¤y;%initialguess=minenergy22tic23xp=l1eqpd(x0,A,[],y,1e¡3);24toc2526norm(f¡xp)/norm(f)27plot(f)28holdon29plot(xp,'r.')30legend('Original','Recovered')TherecoverycanbemadebyusingCVXinsteadofl1-Magic.Justreplacelines19to24by150100150200250300350400450500−3−2−10123OriginalRecoveredFigure1:Script1results.Listing3:UsingCVXfortherecovery.%SolveusingCVX.cvxbeginvariablexp(n);minimize(norm(xp,1));subjecttoA¤xp==y;cvxend1.2UsingaFourierSensingMatrixNowwearegoingtorepeatthesameexperiment,butusingasamplingmatrixbasedontheFourierbasis.Wegeneratethemeasurementmatrixwiththefollowingfunction:Listing4:Fourierbasedmeasurementmatrix.functionA=getAfourier(n,m)tmp=randperm(n);phi=inv(fft(eye(n)));A=phi(tmp(1:m/2),:);A=[real(A);imag(A)];Inordertorecoverthesignalusingl1-magic,nowweneedtousethefunctionl1qc_logbarrierinsteadofl1eq_pd.Ontheotherhand,thereisnoneedtochangeanythingwhensolvingtheproblemwithCVX.SinceingeneralissimplerandclearertouseCVX,weonlyusethisapproachinthefollowingexamples.2SparseSignalintheFrequencyDomainLetstrynowwithasignalsparseinthefrequencydomain.Wegeneratethesignalas:Listing5:Sparsesignalinthefrequencydomain.t=[0:n¡1]';f=cos(2¤pi/256¤t)+cos(2¤pi/128¤t);2LetssolvewitharandomsensingmatrixandCVX.Noticethatnowweneedtospecifytherepresentationbasisª(seeline51):Listing6:RandommeasurementsandCVXrecovery.1A=getArandom(n,m);2y=A¤f;34%SolveusingCVX.5Psi=inv(fft(eye(n)));6cvxbegin7variablexp(n);8minimize(norm(xp,1));9subjectto10A¤Psi¤xp==y;11cvxendFigure2showstheresult.Asexpected,therecoveryisexact.Letsmodifyoursignalslightly,byreplacingoneofthecosinebyasine:0100200300400500600−1.5−1−0.500.511.52OriginalRecoveredFigure2:Recoveryofafrequencydomainsparsesignal.Listing7:Sparsesignalinthefrequencydomain.t=[0:n¡1]';f=cos(2¤pi/256¤t)+sin(2¤pi/128¤t);Figure3showstheresult.Evidentlythereissomethingwrong.TheproblemisthatnowtheFouriercoe±cientshaveanimaginarycomponent,butCVXissearchingforarealx.Thesolutioniseasy,wejustneedtotellCVXtoconsideracomplexx:Listing8:Sparsesignalinthefrequencydomain.%SolveusingCVX.Psi=inv(fft(eye(n)));cvxbeginvariablexp(n)complex;%WeneedtotellCVXthatxpiscomplex!minimize(norm(xp,1));subjecttoA¤Psi¤xp==y;cvxend30100200300400500600−3−2−10123OriginalRecoveredFigure3:Wrongrecovery.Figure4showstheresult.Nowtherecoveryisexact.Wecanalsotrywithahighersparsitylevel.Thefollowingscriptcreateasignalbyaddingsixsinusoidswithrandomperiod,amplitudeandphases.Noticethattheperiodsarechosenfromthevector[163264128256512].Figure5showstheresult.Onceagain,therecoveryisexact.Listing9:Signalmadeof6randomsinuoids.1s=6;2amp=rand(s,1);%amplitudes3periods=[163264128256512];4tmp=randperm(length(periods));5freq=(2¤pi./round(periods(tmp(1:s))));%frequencies6phases=2¤pi¤rand(s,1);7f=zeros(n,1);8t=[0:n¡1]';9fork=1:s,10f=f+amp(k)¤cos(freq(k)¤t+phases(k));11endWenowreplacetheperiodsweareusingby[183264128256512].Noticethattheonlydi®erenceisthatnowthesmallestperiodis18insteadof16.Figure6showstheresult.Nowtheresultisnotexact.Thereasonforthisisthatnowthesignalisnotreallysparse,sinceoneoftheperiodsisnotanintegermultipleofthesignallength.3AcknowledgmentsThankstoDr.MichaelWakinandBorhanSanandajiforhelpingmetosolvesomeoftheissuesIhadwiththecode.1Ifthesignalissparseinthetimedomain,ª=Identitymatrix,thatswhywedidn'tspecifyªinsection1.40100200300400500600−2−1.5−1−0.500.511.52OriginalRecoveredFigure4:SuccessfullyrecoveryaftertellingCVXtouseacomplexvariable.50100150200250300350400450500−3−2−10123OriginalRecoveredFigure5:Recoveryofasignalmadeof6sinusoids.550100150200250300350400450500−3−2−10123OriginalRecoveredFigure6:Recoveryofasignalmadeof6sinusoids.Thesmallestperiodisnow18insteadof16.6