home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 3.ddi / MMLE3.DI$ / ML_DEM1G.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.6 KB  |  87 lines
  1. %ROUTINE TO GEMERATE TESTCASE1 DATA FROM MMLE3 PROGRAM BY MAINE & ILIFFE
  2. %   - PRODUCES   uydata  FOR USE BY   ML_DEMO1 ROUTINE.
  3. %   - STORES DATA IN FILE ML_DEMO1.MAT   SEE LISTING FOR MORE DETAILS
  4.  
  5.  
  6. disp('----------------------------------------------------------------------')
  7. disp('I am about to spend more than fifteen minutes generating 2004 specific')
  8. disp(' random numbers. When done, I will over-write disk file  ml_demo1.mat')
  9. disp(' To stop me enter CTRL-C.  Touch anything else and I go ahead.')
  10. disp('----------------------------------------------------------------------')
  11. pause
  12. disp('starting')
  13.  
  14. uydata=[];
  15. ndp=1001;
  16.  
  17. atrue=-1;btrue=10;ftrue=2;gtrue=1;dt=.01;
  18.  
  19. sqrtq=ftrue*sqrt(dt);% DISCRETE PROCESS NOISE COVARIANCE APPROXIMATION
  20.  
  21. [phi,gam]=c2d(atrue,btrue,dt);
  22.  
  23. t=[0:ndp]'*dt;
  24.  
  25. gauss=[t;t];
  26.  
  27.  
  28. %--------------------------------------------NOISE GENERATION ALGO FROM MMLE3
  29. imod=131071;rdiv=131072;imult=2401;seed=4455;r=[];x=[];y=[];
  30.  
  31. doshift=['seed=round(rem(seed*imult,imod));x=2*seed/rdiv-1;',...
  32.          'seed=round(rem(seed*imult,imod));y=2*seed/rdiv-1;',...
  33.          'r=x^2+y^2;'];
  34.  
  35. for i=1:length(gauss) ;i
  36.      if rem(i,2)==1
  37.         eval(doshift)
  38.            while r>=1
  39.               eval(doshift)
  40.            end
  41.         r=sqrt(-2*log(r)/r);
  42.         gauss(i)=x*r;
  43.      else
  44.         gauss(i)=y*r;
  45.      end
  46. end
  47. % ---------------------------1 dimensional noise vector generation complete
  48.  
  49. statnoise=gauss(1:2:2*ndp-1);   % en
  50. measnoise=[0;gauss(2:2:2*ndp)]; % (Eta)
  51. clear gauss
  52.  
  53.  
  54. u=rem(fix(t+.005),2);
  55.  
  56. %Maine and Iliff's model is
  57.  
  58. %        X(n+1)=PHI * X(n) + GAM * [U(n)+U(n+1)]/2
  59.  
  60. %  This average value hold avoids the half-sample delay in the usual ZOH.
  61. %  We must thus average the input samples to reproduce their data.
  62.  
  63.    u_data(1:ndp,1)=(u(1:ndp)+u(2:ndp+1))/2;  % Average the inputs (see end)
  64.  
  65. xnonoise=dlsim(phi,gam,1,0,u_data); %  Simulation with no noise
  66.  
  67. xnosig=dlsim(phi,1,1,0,statnoise); %  Separately simulate the process noise
  68.  
  69. xnosig=xnosig(1:ndp);measnoise=measnoise(1:ndp);
  70.  
  71. y_data = xnonoise + xnosig*sqrtq + measnoise*gtrue; %  Combine the components
  72.  
  73. y_data=y_data(1:ndp);u_data=u_data(1:ndp);
  74.  
  75. uydata=[u_data y_data];
  76.  
  77. save ml_demo1   uydata   xnonoise   xnosig   measnoise    statnoise
  78.  
  79. clear
  80.  
  81. % The output data are uydata=[u_data y_data]. The mmle routine in PC-Matlab
  82. % doesn't internally average inputs as Maine & Iliff's MMLE3 program does,
  83. % so u_data is the correct input to use. MMLE3 would use u.
  84.  
  85. disp('--------------------------------ml_dem1g.m finished.   Stack cleared.')
  86.  
  87.