ls.m

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

function [t_array,f_array,a_array,phi_array,select_array,residual] = LS(t,s,t_scan_num,f_scan_num,threshold,amp,f_width,delt,k,ebs)
    % this function is to get the least-square wavelets
    %t sample time
    %s sample signal
    %t_scan_num   better to be odd number
    %f_scan_num   better to be odd number
    %threshold    control when to stop decomposite
    %amp          etermine the amplitude of the morlet window  default:=1
    %f_width      determine the wavenumber of the window default: =0.09
    %delt         determine the breath of the window default: =0.5*1.0e-19
    %k            constant value that controls the wavelet breadth  
    %ebs          a small number which makes the solution stable. default:0
    i=sqrt(-1);   
    [m,n]=size(t);  % n: sample number
    Ts=(t(n)-t(1))/(n-1);
    fs=1/Ts;
    s=hilbert(s);% Hilbert Transform
    residual=s;
    E=sum((abs(s)).^2);
    E_threshold=threshold*E;  
    atom_num=0;
    E_gap=E;
    while E_gap>E_threshold
        atom_num=atom_num+1;
        %---------------- find center time & peak frequency &&improve accuracy
        [amplitude_max,j]=max(abs(residual));
        t0_center=t(j); 
        if t_scan_num==1
            t_scan=t0_center;
        else
            t_scan_range=Ts;  %automatically define
            t_scan=linspace((t0_center-t_scan_range/2),(t0_center+t_scan_range/2),t_scan_num);
        end
        E_min=E_gap;  %residual energy minimum
        %time scan
        for ii=1:t_scan_num
            t_center=t_scan(ii);
            [window,x]=morlet(t(1),t(n),n,t_center,delt,amp,f_width);
            residual_short=residual.*window;
            NFFT = 2^nextpow2(n);
            spectrum=fft(residual_short,NFFT)/n;
            spectrum_abs=2*abs(spectrum);
            f = fs*linspace(0,1,NFFT);
            [spectrum_abs_max,j]=max(spectrum_abs(1:NFFT));
            f0_peak=f(j);
            if f_scan_num==1
                f_scan=f0_peak;
            else
                f_scan_range=fs/NFFT;  %///////////////how to automatically define?
                f_scan=linspace((f0_peak-f_scan_range/2),(f0_peak+f_scan_range/2),f_scan_num);
            end
            %frequency scan 
            for jj=1:(f_scan_num)
                f_peak=f_scan(jj);
                %---------------find amplitude &phase
                f_i=f_peak;
                t0=t_center;
                [psi_real,x]=morlet_cos_yueqin_st(t(1),t(n),n,f_i,t0,k);
                [psi_imag,x]=morlet_sin_yueqin_st(t(1),t(n),n,f_i,t0,k); 
                wavelet=(psi_real)'+i*(psi_imag)';
                [n1,n2]=size(wavelet);
                %---------------compute A
                wavelet_H=wavelet';
                wavelet_E=wavelet_H*wavelet;
                I=eye(n2);
                A=(inv(wavelet_E+ebs*I))*wavelet_H*residual.';
                a=abs(A);
                phi=angle(A);
                %judge if to break
                select=A*wavelet.';
                residual_temp=residual-select;
                E=sum((abs(residual_temp)).^2);
                if E<E_min
                    E_min=E;
                    t_array(atom_num)=t_center;
                    f_array(atom_num)=f_peak;
                    a_array(atom_num)=a;
                    phi_array(atom_num)=phi;
                    select_array(atom_num,1:n)=select;
                    residual_next=residual_temp;
                end
            end
        end
        residual=residual_next;
        E_gap=E_min;
        %title('2 Iteration')
        %xlabel('Time(s)')
        %ylabel('Amplitude')
        %axis([time(1) time(nb) -200 300]);
        %figure_num=figure_num+1;
        %figure (figure_num)
    end