利用FFT实现MP的稀疏分解参考程序

附录D利用FFT实现基于MP的信号稀疏分解参考程序%%MATLABscriptforMPalgorithmbasedonFFT%%%**********************Preparationpart**********************%%**********************part1:imputsignal**********************%%readsignaldata(tobeprocessed)%Herewecaninputdifferenttypesofdata,tochecktheprogram.%Inthisexample,wereadonesignaldatafromdisk,whichissavedinthe%format:savebat.datbat-asciiloadbat.dat_ascii;signal=bat;%**********************part2:setparameters*********************************iterative_number=30;%thematchingpursuitprocessingiterativenumber[a,N]=size(signal);%thelengthofsignalandthelengthofatoms;Nsignal_reconstruct=zeros(1,N);signal_r=signal;a_base=2;%scale:thescaleisdeterminedbyaandj;j_min=0;j_max=log2(N);u_base=l/2;%thetransmissionordisplacementushoulddeterminedasfollowingv_base=pi;%thefrequencyvk_min=0;w_base=pi/6;%thephasewi_min=0;i_max=12;%*******************wipeoffthedirectcurrentvector*****************signal_reconstruct=(1/N)*sum(signal);signal_r=signal-signal_reconstruct;%*************createdictionary*************************************[atoms]=dic_a(signal_r,N,a_base,j_min,j_max,u_base,v_base,k_min,w_base,i_min,i_max%*************startcalculation*************************************forn=1:iterative_numbern%*******thefollowingprogramusesonesubroutinetoselectthebestatom*********[proj,scale,translation,freq,phase]=select_best(signal_r,N,a_base,j_min,j_max,u_base,v_base,k_min,w_base,i_min,i_max,atoms);%*****reconstructthebestatomfromtheparametersgottedbyabovesubroutine********t=0:N-1;t=(t-translation)/scale;g1=(1/sqrt(scale))*exp(-pi*t.*t).*cos(freq*t+phase);g=g1/sqrt(sum(g1.*g1));%************reconstructsignalandrenewtheresidual*******************signal_reconstruct=signal_reconstruct+proj*g;signal_r=signal_r-proj*g;%************thefigureofresult*******************subplot(221);plot(signal);%originalsiganlsubplot(222);plot(g);%theatomsubplot(223);plot(signal_r);%theresidualsubplot(224);plot(signal_reconstruct);%reconstructsignalend%************************endofmainprogramfile**************%*****thefollowingistheatomdictionarycreationprogram******function[atoms]=dic_a(signal_r,N,a_base,j_min,j_max,u_base,v_base,k_min,w_base,i_min,i_max)%thissubroutineistocreatethedictionary%INPUT%thesignal_r:thesignalortheresidualofthesignaltobedecomposed%theN:thelengthofthesignaloroftheresidualofthesignalorthelengthoftheatoms%parameters:theparametertoconstructthedictionary,ithasmuchinfluenceonthe%speedofthedecomposition%OUTPUT%atoms:thedictionarysize_dic=1;atoms=zeros(size_dic,N);%************************creatdictionary**************forj=j_min:j_maxfork=k_min:2^(j+1)fori=i_min:i_maxsize_dic=size_dic+1;%themunberofatomss=a_base^j;%scaleoftheatomu=N/2;%transmissionoftheatomv=k*(1/s)*v_base;%frequencyoftheatomw=i*w_base;%phaseoftheatomt=0:N-1;t=(t-u)/s;g=(1/sqrt(s))*exp(-pi*t.*t).*cos(v*t+w);%generatetheatomg=g/sqrt(sum(g.*g));%normalizationatoms(size_dic,:)=g;endendend%***************endofthedictionarycreationsubroutineprogram********%*****thefollowingisthebestatomselectingsubroutineprogramusingFFT******Function[projscale,translation,freq,phase]=select_best(signal_r,N,a_base.j_minj.u^base,v_base.k_min,w_base,i_min,i_max,atoms)%thissubroutineistoselectinthedictionarythebestatomsuitedthesiganlorthe%residualofthesignalusingFFTtoreducecomputationaltime%INPUT%thesignal_r:thesignalortheresidualofthesignaltobedecomposed%theN:thelongthofthesignaloroftheresidualofthesignalorthelengthoftheatoms%parameters:theparametertoconstructthedictionary,ithasmuchinfluenceonthe%speedofthedecomposition%OUTPUT%proj:theprojectionofthesignalortheresidualofthesignalonthebestatom%thescale:thescaleofthebestatom(sintheformula)%thetranslation:thetranslationofthebestatom(uintheformula)%thefreq:thefrequencyofthebestatom(vintheformula)%phase:thephaseofthebestatom(wintheformula)%proj_trans:todeterminewhichprojectionisbiggest%************************intializtion***********************[mn]=size(atoms);re=zeros(5,m);proj_trans=0;proj=0;size_dic=0;%************************selcetthebestatom***********************forj=j_min:j_maxfork=k_min:2^(j+1)fori=i_min:12size_dic=size_dic+1;g=[atoms(size_dic,N:-1:1)];co=real(ifft(fft(signal_r,2*N).*fft(g,2*N)));recon=co((N/2+1):3*N/2);[proj_transe]=max(abs(recon));proj_trans=rccon(e);%choosethebiggestprojectionre(:,size_dic)=[proj_transjeki];endendendse=re(1,:);[proj_transe]=max(abs(se));proj=se(e);scale=re(2,e);translation=re(3,e);freq=re(4,e);%*************computetheparameteroftheselectedatom********************phase=re(5,e);scale=a_base^scale;freq=freq*(1/scale)*v_base;phase=phase*w_base;