xgmp2.m

GPS software toolbox for GPS receiver development

%                                xgmp2.m
%  Scope:   This MATLAB program generates second order Gauss-Markov sequence, 
%           plots the generated sequence and determines and plots the 
%           normalized auto-correlation sequence.
%  Usage:   xgmp2
%  Inputs:  - number of steps in the sequence to be generated (nstep)
%           - the value of the natural frequency in radians/second
%           - the value of the damping factor (less than 1.)
%           - the value of the constant c**2 in meters**2 
%           - the value of the time step in seconds
%           - value of the initial seed, default value is 0
%           - name of the output file, if the generated data are saved
%  Outputs: - plot of the second order Gauss-Markov sequence - first state
%             component
%           - plot of the second order Gauss-Markov sequence - second state
%             component
%           - plot of the normalized auto-correlation sequence - first state
%             component (optional)
%           - plot of the normalized auto-correlation sequence - second state
%             component (optional)
%  Remark:  A default set of input data is provided.
%  External Matlab macros used: genrn, gmp2, rms, xcorr (from Signal 
%                               Processing Toolbox) - optional
%  Last update:  06/28/00
%  Copyright (C) 1996-00 by LL Consulting. All Rights Reserved.

clear
close all
yes = 'y';

%  Initialization - input data for second order Gauss-Markov process

disp('  ');
disp('Specify the number of steps in the sequence.');
nstep = input('Example : 256, 512, 1024, 2048, 4096;  Make selection --> ');
disp(' ');
answer1 = input('Do you want to use the default data? (y/n)[y] ','s');
if  isempty(answer1)
   answer1 = yes;
end
if (strcmp(answer1,yes) == 1)
   w0 = 0.012;                %  natural frequency in radians/second
   beta = 0.707106781;        %  value of  1/sqrt(2.0) , damping factor
   csq = 0.002585;            %  constant  c**2  in meters**2          
   deltat = 1.0;              %  time step in seconds
else
   w0 = input('Specify the natural frequency in radians/second -->  ');
   beta = input('Specify the damping factor (less than 1.) -->  ');
   csq = input('Specify the constant c**2 in meters**2 -->  ');
   deltat = input('Specify the time step in seconds -->  ');
end
disp('  ');
disp('Select the initial seed value, the default is 0');
answer2 = input('Do you want to use the default value? (y/n)[y] ','s');
if  isempty(answer2)
   answer2 = yes;
end
if (strcmp(answer2,yes) == 1)
   iseed = 0;
else
   iseed = input('Specify the value --> ');
end  
rand('seed',iseed);
disp('  ');

%  Generate the second order Gauss-Markov sequence

[xp,xv] = gmp2(nstep,w0,beta,csq,deltat);

%  Determine mean, standard deviation and root mean square (rms) for the 
%  generated second order Gauss-Markov sequence - first state component

xmean = mean(xp);
xstd = std(xp);
xrms = rms(xp);
temp1 = ['mean = ',num2str(xmean)];
temp2 = ['st.dev. = ',num2str(xstd)];
temp3 = ['rms = ',num2str(xrms)];
g = [temp1,' ;   ',temp2,' ;   ',temp3];

%  Plot the second order Gauss-Markov time sequence - first state component

fff = num2str(nstep);
dd = 'Second order Gauss-Markov process - first state component ';
ee = 'First state component';
ff = 'Time sequence (Total sequence length = ';
ff = [ff fff ')']; 
disp('Execute second order Gauss-Markov process graph - first component -');
disp('Select the mouse position to insert text on the graph.');
disp('Press a key to continue...');
pause
figure(1)
subplot(2,1,1)
plot(xp),...
title(dd), ylabel(ee), xlabel(ff), grid,...
gtext(g);                   %  mouse placement of the text on a graph
disp('  ');
disp('Executed graph - check the window associated to this figure.');
disp('Press a key to continue...');
pause

%  Determine mean, standard deviation and root mean square (rms) for the 
%  generated second order Gauss-Markov sequence - second state component

xmean = mean(xv);
xstd = std(xv);
xrms = rms(xv);
temp1 = ['mean = ',num2str(xmean)];
temp2 = ['st.dev. = ',num2str(xstd)];
temp3 = ['rms = ',num2str(xrms)];
g = [temp1,' ;   ',temp2,' ;   ',temp3];

%  Plot the second order Gauss-Markov time sequence - second state component

dd = 'Second order Gauss-Markov process - second state component';
ee = 'Second state component';
disp('  ');
disp('Execute second order Gauss-Markov process graph - second component -');
disp('Select the mouse position to insert text on the graph.');
disp('Press a key to continue...');
pause
figure(1)
subplot(2,1,2)
plot(xv),...
title(dd), ylabel(ee), xlabel(ff), grid,...
gtext(g);                   %  mouse placement of the text on a graph
disp('  ');
disp('Executed graph - check the window associated to this figure.');
disp('Press a key to continue...');
pause

%  Plot the normalized auto-correlation sequence

disp('  ');
disp('For the auto-correlation sequence (first/second state component) -  ');
disp('   the function XCORR from Signal Processing Toolbox is required.');
answer3 = input('Plot the auto-correlation sequence? (y/n)[y] --> ','s');
if  isempty(answer3)
   answer3 = yes;
end  
disp('  ');

%  Plot the normalized auto-correlation sequence - first state component

if strcmp(answer3,yes) == 1
   b = -(nstep-1):1:(nstep-1);
   a = xcorr(xp','coeff');
   aa = 'Second order Gauss-Markov process - first state component';
   bb = 'Normalized auto-correlation sequence';
   cc = 'Sequence lag  (Total sequence length = ';
   cc = [cc fff ')'];
   figure(2)
   subplot(2,1,1)
   plot(b,a'),...
   title(aa), ylabel(bb), xlabel(cc), grid;
end   

%  Plot the normalized auto-correlation sequence - second state component

if strcmp(answer3,yes) == 1
   a = xcorr(xv','coeff');
   aa = 'Second order Gauss-Markov process - second state component';
   bb = 'Normalized auto-correlation sequence';
   figure(2)
   subplot(2,1,2)
   plot(b,a'),...
   title(aa), ylabel(bb), xlabel(cc), grid;
   disp('Auto-correlation graph - check the window associated to this figure.');
   disp('Press a key to continue...');
   disp('  ');
   pause
end

%  Save the generated data into a specified file if desired

answer4 = input('Do you want to save the generated data? (y/n)[y] --> ','s');
if  isempty(answer4)
   answer4 = yes;
end
if strcmp(answer4,yes) == 1
   disp('  ');
   f2 = input('Specify the output filename  -->  ','s');
   for  k = 1:nstep
      fprintf(f2,'%24.14e    %24.14e\n',xp(k),xv(k));
   end
end

disp('  ');
disp('End of the program XGMP2');
disp('  ');