smoothing.m

利用gps多天线载波相位技术

%%========================================
%%     Toolbox for attitude determination
%%     Zhen Dai
%%     [e-mail]
%%     ZESS, University of Siegen, Germany
%%     Last Modified  : 1.Sep.2008
%%========================================
%% Functions:
%%      Smooth the code data using carrier phase data
%% Input parameters:
%%      mRawCode-> Raw code data obtained directly from RINEX file
%%      mRawPhase-> Raw phase data obtained directly from RINEX file
%%      smooth_interval -> Smoothing interval in epochs 
%% Output:
%%      mCodeSmoothed ->Smoothed code data 
%% Remarks:
%%      1. Applied for single frequency data
%%      2. No strict cycle-slip detection and correction. Only a rough
%%      check for cycle-slips. Pre-processed phase data are therefore
%%      expected. 
%% References:
%%      1. R.Hatch: "The Synergism of GPS Code and Carrier Measurements"
%%          International Geodetic Symposium on Satellite Doppler
%%          Positioning, 3rd, Las Cruces, NM, February 8-12, 1982, vol. 2 , pp. 1213-1231.
%%      2. B.Hofmann-Wellenhof, H.Lichtenegger and J.Collins: GPS Theory
%%          and practice. 2001. Fifth revised edition. Springer, Wien, New York.
%%          pp. 94-97.  

function mCodeSmoothed=Smoothing(mRawCode,mRawPhase,smooth_interval)
 
lambda1=0.19029;
 
waitbar_smooth = waitbar(0,'Smoothing code data');

totalepoch=size(mRawCode,2);
totalgps=32; %% assuming there are totally 32 GPS satellites

for satid=1:1:totalgps,
    for epoch=1:1:totalepoch, 
        %% current processed epoch
        current_epoch=epoch;
        %% starting epoch for smoothing, for example, 
        %% if  smoothing interval is 100 epochs,only the data of 100 
        %% epochs priori to the current epoch are considered.
        first_smooth_interval=current_epoch-smooth_interval+1;
        %% for the case of the beginning of data
        if first_smooth_interval<=0, 
            first_smooth_interval=1;
        end
        %% vR0 is a temp vector containing R(t1)i, namely the  temporary
        %% smoothed code of the epochs between the starting epoch
        %% and the current epoch
        clear vR0
        %% If no phase avaliable, then the smoothing procedure restarts.
        %% Otherwise, accumulate the data
        if mRawPhase(satid, current_epoch)==0, %% no phase data
            first_smooth_interval=current_epoch+1;
            mCodeSmoothed(satid, current_epoch)=mRawCode(satid, current_epoch);
        else %% Phase data avaliable
            %% Accumulate the data for furture average.
            t=0;
            for i=first_smooth_interval:1:current_epoch,
                t=t+1; 
                vR0(t)=mRawCode(satid,i)-mRawPhase(satid, i)*lambda1;
            end
            %% Average the previous data and updated with the current ones
            mCodeSmoothed(satid, current_epoch)=mean(vR0)+(mRawPhase(satid, current_epoch))*lambda1;
            %% if the difference between smoothed and unsmoothed data are larger than 5 meters
            %% it might probably be a cycle-slip, so reset the smoothing.
            if abs(mCodeSmoothed(satid, current_epoch)-mRawCode(satid, current_epoch))>5,
                mCodeSmoothed(satid, current_epoch)=mRawCode(satid, current_epoch);
                first_smooth_interval=current_epoch+1;
            end
        end        
    end %% Loop for epochs
    str=sprintf('Code data of satellite PRN. %d have been smoothed',satid);
    waitbar(satid/totalgps,waitbar_smooth,str)
end %% Loop for satellite
close(waitbar_smooth);