mp_least_square_synthetic.m

基于morlet小波的最小二乘匹配追踪算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%mp-least square algorithm
% 1 definiton
clear all
i=sqrt(-1);
lb=-0.4*10^(-6); ub=0.4*10^(-6); nb=1024;  %nb should be 10^n;
Ts=((ub-lb)/(nb-1));
fs=1/Ts;
figure_num=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2 create synthetic signal 
out2=(0:nb-1)*Ts+lb;   %out2 = linspace(lb,ub,nb);
f=200*10^6;
a=0.2;
phi=pi/4;
center_time=0.00*10^(-6);
atom1 = a*exp(-((out2-center_time).^2)*(f^2)*log(2)) .* cos(2*pi*f*(out2-center_time)+phi)+i*a*exp(-((out2-center_time).^2)*(f^2)*log(2)) .* sin(2*pi*f*(out2-center_time)+phi);
f=40*10^6;
a=0.7;
phi=pi/3;
center_time=0.20*10^(-6);
atom2 = a*exp(-((out2-center_time).^2)*(f^2)*log(2)) .* cos(2*pi*f*(out2-center_time)+phi)+i*a*exp(-((out2-center_time).^2)*(f^2)*log(2)) .* sin(2*pi*f*(out2-center_time)+phi);
f=80*10^6;
a=0.5;
phi=pi/2;
center_time=0.3*10^(-6);
atom3 = a*exp(-((out2-center_time).^2)*(f^2)*log(2)) .* cos(2*pi*f*(out2-center_time)+phi)+i*a*exp(-((out2-center_time).^2)*(f^2)*log(2)) .* sin(2*pi*f*(out2-center_time)+phi);
f=60*10^6;
a=0.9;
phi=0;
center_time=-0.2*10^(-6);
atom4 = a*exp(-((out2-center_time).^2)*(f^2)*log(2)) .* cos(2*pi*f*(out2-center_time)+phi)+i*a*exp(-((out2-center_time).^2)*(f^2)*log(2)) .* sin(2*pi*f*(out2-center_time)+phi);
signal_original=atom1+atom2+atom3+atom4;
figure (figure_num)
plot(out2,real(signal_original),'k-');
%figure_num=figure_num+1;
%figure (figure_num)
hold on
%plot(out2,imag(signal_original),'r-');
C1=real(signal_original);
D1=imag(signal_original);
%signal=signal_original;
%plot(out2,abs(signal_original),'r-');
%%%%%%%%%%%%%%%%%%%%%%%%add guassion noise
[m,n]=size(signal_original);
%noise=1.0e-2*randn(m,n); 
%noise=hilbert(noise);
%signal=signal_original+noise;
SNR=-1;
signal = awgn(signal_original,SNR,'measured');   %add noise
figure_num=figure_num+1;
figure (figure_num)
plot(out2,real(signal),'r-');
%hold on
%plot(out2,imag(signal),'r-');

%plot(out2,abs(signal),'g-');
%plot(out2,imag(signal),'r-');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ST time-frequency distribution
%[st,t,f] = st(real(signal),0,3*1.0e6,Ts,1);
%st=fftshift(st,2);
%st_abs=abs(st);
%figure_num=figure_num+1;
%figure (figure_num)
%contourf(out2,f,st_abs);
%grid on
%xlabel('t(s)');
%ylabel('f(Hz)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 mp circulation
%signal=signal_original;
residual=signal;
energy=sum((abs(signal)).^2);
energy_threshold=0.005*energy;  %///////////////how to automatically define?
atom_num=0;
energy_gap=energy;
while energy_gap>energy_threshold
    atom_num=atom_num+1;
    % 4.1 find center time & peak frequency &&improve accuracy
    [amplitude_max,j]=max(abs(residual));
    %max=find(diff(sign(diff(abs(residual))))== - 2)+1;
    center_time_0=out2(j);           
    t_scan_range=Ts;  %automatically define
    t_scan_num=9;     %///////////////how to automatically define?
    t_scan=linspace((center_time_0-t_scan_range/2),(center_time_0+t_scan_range/2),t_scan_num);
    energy_min=energy_gap;  %residual energy minimum
    for ii=0:(t_scan_num-1)
        center_time=center_time_0-t_scan_range/2+ii*t_scan_range/(t_scan_num-1);
        delt=1.0e-18; %determine the breath of the window 
        amp=1;
        f_width=0.09;
        [psi,x]=morlet(lb,ub,nb,center_time,delt,amp,f_width);
        %figure_num=figure_num+1;
        %figure (figure_num)
        %plot(x,psi,'m-');
        hold on
        residual_short=residual.*psi;
        %figure_num=figure_num+1;
        %figure (figure_num)
        %plot(out2,residual_short);
        %figure_num=figure_num+1;
        NFFT = 2^nextpow2(nb);
        spectrum=fft(residual_short,NFFT)/nb;
        spectrum_abs=2*abs(spectrum);
        f = fs*linspace(0,1,NFFT);
        %figure (figure_num)
        %plot(f,2*spectrum_abs(1:NFFT));
        %figure_num=figure_num+1;
        [spectrum_abs_max,j]=max(spectrum_abs(1:NFFT));
        frequency_peak_0=f(j);
        %f_scan_range=fs/4*10^(-2);  %///////////////how to automatically define?
        f_scan_range=fs/NFFT;  %///////////////how to automatically define?
        f_scan_num=3;                   %///////////////how to automatically define?
        f_scan=linspace((frequency_peak_0-f_scan_range/2),(frequency_peak_0+f_scan_range/2),f_scan_num);
        %should add frequency scan here later!
        k=1;  %constant value that controls the wavelet breadth  
        for jj=0:(f_scan_num-1)
            frequency_peak=frequency_peak_0-f_scan_range/2+jj*f_scan_range/(f_scan_num-1);
            % 4.2 find amplitude &phase
            f_i=frequency_peak;
            t0=center_time;
            [psi_real,x]=morlet_cos_yueqin_st(lb,ub,nb,f_i,t0,k);
            if atom_num==2;
                %plot(x,psi_real,'r-');
            end
            %figure (figure_num)
            %plot(x,psi_real);
            %figure_num=figure_num+1;
            [psi_image,x]=morlet_sin_yueqin_st(lb,ub,nb,f_i,t0,k); 
            %figure (figure_num)
            %plot(x,psi_image);
            %figure_num=figure_num+1;
            wavelet_1=psi_real+i*psi_image;
            wavelet_1=wavelet_1.';
            [n,m]=size(wavelet_1);
              %compute A
            ebs=0; %a small number which makes the solution stable 
            wavelet_1_H=wavelet_1';
            %figure_num=figure_num+1;
            %figure (figure_num)
            %plot(out2,real(wavelet_1_H),'k-');
            hold on
            %plot(out2,imag(wavelet_1_H),'r-');
           
            wavelet_temp=wavelet_1_H*wavelet_1;
            I=eye(m);
            temp1=wavelet_temp+ebs*I;
            temp2=inv(temp1);
            
            %plot(out2,A_2_real_temp-A_1_real_temp,'B-');
           % delt_temp=sum(A_2_real_temp-A_1_real_temp);
            hold on
            %plot([out2(1),delt_temp],[out2(nb),delt_temp],'r-');
           
            A=temp2*wavelet_1_H*residual.';
            a=abs(A);
            phi=angle(A);
            %judge if to break
            substraction_temp=A*wavelet_1.';
            residual_temp=residual-substraction_temp;
            energy=sum((abs(residual_temp)).^2);
            if energy<energy_min
                energy_min=energy;
                center_array(atom_num)=center_time;
                frequency_array(atom_num)=frequency_peak;
                a_array(atom_num)=a;
                phi_array(atom_num)=phi;
                substraction(atom_num,1:nb)=substraction_temp;
                residual_next=residual_temp;
            end
        end
    end
    residual=residual_next;
    plot(out2,substraction(atom_num,1:nb),'r-');
    hold on
   plot(out2,residual,'b-');
    axis([out2(1) out2(nb) -0.8 1]);
    energy_gap=energy_min;
end
%compute error & draw the restruction
restruction=zeros(1,nb);
for ii=1:atom_num
    restruction=restruction+substraction(ii,1:nb);
end
error=signal-restruction;
plot(out2,real(restruction),'b-');
hold on
plot(out2,real(error),'r-');
figure_num=figure_num+1;
figure (figure_num)
error_dB=20*log10(real(error./signal));
plot(out2,error_dB,'b-');
% delete below
figure_num=figure_num+1;
figure (figure_num)
plot(1:atom_num,center_array,'k*')