mp_least_square_gpr.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%mp-least square algorithm
% 1 definiton
clear all
i=sqrt(-1);
%import GPR signal
receiver_1=load('e:\gpr_test\case 1\receiver_1.txt');
out2=receiver_1(2:(receiver_1(1,1)+1),1);
signal_real=receiver_1(2:(receiver_1(1,2)+1),2);
[nb,a]=size(out2);
Ts=(out2(nb)-out2(1))/(nb-1);
fs=1/Ts;
figure_num=1;
figure (figure_num)
%plot(out2,signal_real,'k-');
hold on
% Hilbert Transform
signal=hilbert(signal_real);
%plot(out2,real(signal),'r-');
grid on
signal=signal.';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%WVD distribution
%[tfr,t1,f1] = tfrwv(signal_real);
%tfr=fftshift(tfr,1);
%f=fs*linspace(0,1,nb);
%figure_num=figure_num+1;
%figure (figure_num)
%contourf(t,f,tfr);
%grid on
%xlabel('t(s)');
%ylabel('f(Hz)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 mp circulation
residual=signal;
energy=sum((abs(signal)).^2);
energy_threshold=0.01*energy;  %///////////////how to automatically define?
atom_num=0;
energy_gap=energy;
figure_num=figure_num+1;
figure (figure_num)
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));
    %amplitude_max_id=find(diff(sign(diff(real(residual))))== - 2)+1;
    center_time_0=out2(j);     
    %plot([out2(j),out2(j)],[0,300],'k-');
    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);
    %t_scan(t_scan_num+1)=center_time_0;  % include the original center_time 
    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=0.5*1.0e-19; %determine the breath of the window 
        amp=1;  %determine the amplitude of the window
        f_width=0.09;
        [psi,x]=morlet(out2(1),out2(nb),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,'g-');
        %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);
        %spectrum_abs_max_id=find(diff(sign(diff(abs(residual))))== - 2)+1;
        %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=0.02;  %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(out2(1),out2(nb),nb,f_i,t0,k);
            if atom_num==2;
                %plot(x,psi_real,'r-');
            end
            %figure_num=figure_num+1;
            %figure (figure_num)
            %plot(x,psi_real);
            %figure_num=figure_num+1;
            [psi_image,x]=morlet_sin_yueqin_st(out2(1),out2(nb),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';
            wavelet_temp=wavelet_1_H*wavelet_1;
            I=eye(m);
            temp1=wavelet_temp+ebs*I;
            temp2=inv(temp1);
            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-');
    energy_gap=energy_min;
    %title('2 Iteration')
   % xlabel('Time(s)')
    %ylabel('Amplitude')
    axis([out2(1) out2(nb) -200 300]);
    %figure_num=figure_num+1;
    %figure (figure_num)
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;
figure (figure_num)
plot(out2,real(restruction),'b-');
hold on
plot(out2,error,'r-');
figure_num=figure_num+1;
figure (figure_num)
error_dB=20*log10((real(error)./real(signal)));
plot(out2,error_dB,'b-');
% delete below
figure_num=figure_num+1;
figure (figure_num)
plot(1:atom_num,center_array,'k*')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
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.01*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.2*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=atom1+atom2+atom3+atom4;
figure (figure_num)
plot(out2,real(signal),'k-');
%figure_num=figure_num+1;
hold on
%plot(out2,imag(signal),'r-');