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

附录D利用FFT实现基于MP的信号稀疏分解参考程序% MATLAB script for MP algorithm based on FFT%* Preparation part *%* part 1:imput signal *%read signal data (to be processed)% Here we can input different types of data, to check the program.% In this example, we read one signal data from disk, which is saved in the% format: save bat.dat bat -asciiload bat.dat _ascii;signal=bat;%*part 2: set parameters *iterative_number= 30 ; % the matching pursuit processing iterative numbera,N=size(signal); % the length of signal and the length of atoms; Nsignal_reconstruct=zeros(1,N);signal_r=signal;a_base=2; % scale: the scale is determined by a and j;j_min=0;j_max=log2(N); u_base=l/2; % the transmission or displacement u should determined as followingv_base=pi; % the frequency vk_min= 0 ;w _base=pi/ 6; % the phase wi_min= 0 ;i_max= 12 ;%*wipe off the direct current vector *signal_reconstruct=( 1 /N)*sum(signal);signal_r=signal-signal_reconstruct;%* create dictionary *atoms=dic_a(signal_r,N,a_base,j_min,j_max,u_base,v_base,k_min,w_base,i_min,i_max%* start calculation *for n=1: iterative_numbern% *the following program uses one subroutine to select the best atom*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);%*reconstruct the best atom from the parameters gotted by above subroutine*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);%* reconstruct signal and renew the residual*signal_reconstruct=signal_reconstruct+proj*g;signal_r=signal_r-proj*g;%*the figure of result *subplot(221);plot(signal); % original siganlsubplot( 222 );plot(g); %the atomsubplot( 223 );plot(signal_r); %the residualsubplot( 224 );plot(signal_reconstruct); %reconstruct signalend% * end of main program file *% * the following is the atom dictionary creation program*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)% this subroutine is to create the dictionary% INPUT% the signal_r:the signal or the residual of the signal to be decomposed% the N: the length of the signal or of the residual of the signal or the length of the atoms% parameters :the parameter to construct the dictionary, it has much influence on the%speed of the decomposition%OUTPUT% atoms:the dictionarysize_dic=1;atoms=zeros(size_dic,N);% * creat dictionary *for j=j_min:j_maxfor k=k_min:2(j+1)for i=i_min:i_maxsize_dic=size_dic+1; %the munber of atomss=a_basej; %scale of the atomu=N/2; % transmission of the atomv=k*(1/s) * v_base; %frequency of the atomw=i*w_base; %phase of the atomt=0:N-1;t=(t-u)/s;g=(1/sqrt(s)*exp(-pi*t.*t).*cos(v*t+w); %generate the atomg=g/sqrt(sum(g. *g); %normalizationatoms(size_dic,:)=g;endendend%* end of the dictionary creation subroutine program*% *the following is the best atom selecting subroutine program using FFT*Function projscale,translation,freq,phase=select_best(signal_r,N,a_base.j_minj.ubase, v_base.k_min,w_base,i_min,i_max,atoms)% this subroutine is to select in the dictionary the best atom suited the siganl or the %residualof the signal using FFT to reduce computational time% INPUT% the signal_r:the signal or the residual of the signal to be decomposed% the N: the longth of the signal or of the residual of the signal or the length of the atoms% parameters :the parameter to construct the dictionary , it has much influence on the%speed of the decomposition%OUTPUT% proj: the projection of the signal or the residual of the signal on the best atom% the scale: the scale of the best atom (s in the formula)% the translation : the translation of the best atom (u in the formula)% the freq:the frequency of the best atom (v in the formula)% phase: the phase of the best atom (w in the formula)% proj_trans :to determine which projection is biggest%*intializtion* m n=size(atoms);re=zeros(5,m);proj_trans=0;proj=0;size_dic=0;%*selcet the best atom*for j=j_min:j_maxfor k=k_min:2(j+1)for i=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_trans e=max(abs(recon);proj_trans=rccon(e); %choose the biggest projectionre(:,size_dic)=proj_trans j e k i;endendendse=r