fig_16_3_16_4_example_demo_profile

% Fig. 16.3 page 731 and 16.4 page 732
 demo profile reconstruction with differente priors.

 Wolfgang Förstner 2014-10-05
 last changes: Susanne Wenzel 09/16
 wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de

0001 %% Fig. 16.3 page 731 and 16.4 page 732
0002 % demo profile reconstruction with differente priors.
0003 %
0004 % Wolfgang Förstner 2014-10-05
0005 % last changes: Susanne Wenzel 09/16
0006 % wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de
0007 
0008 addpath(genpath('../../General-Functions/'));
0009 addpath('Functions');
0010 
0011 close all
0012 
0013 %% plot settings
0014 ss = plot_init;
0015 
0016 %% initialize random number generation by fixed seed
0017 init_rand_seed(4);
0018 
0019 disp(' --------------------------------------------------')
0020 disp(' ------ Fig. 16.3 demo profile reconstruction -----')
0021 disp(' --------------------------------------------------')
0022 
0023 
0024 %% generate true profile
0025 
0026 disp(' ------ generate true profile -----')
0027 N = 200;   % number of points total
0028 N1 = 70;   % number of points in first (flat) part of the profil
0029 N2 = N-N1; % remaing points
0030 
0031 % first AR(1) - flat part of the profil
0032 sigma_e = 1.0;  % process noise
0033 sigma_n = 0.5;  % observation noise
0034 a = 0.98;       % decay
0035 
0036 x1 = zeros(N1,1);   % true profile
0037 y1 = zeros(N1,1);   % observed sample points
0038 for n=1+1:N1
0039     x1(n) = sum(a.*x1(n-1:-1:n-1))+randn(1)*sigma_e;
0040     % discontinuity at position 50
0041     if n==50
0042         x1(n)=x1(n)+8;
0043     end
0044     y1(n) = x1(n)+randn(1)*sigma_n;
0045 end
0046 % enforce slope 0
0047 for i=1:N1
0048     y1(i)=y1(i) -  (i-1)*(x1(N1)-x1(1))/(N1-1);
0049     x1(i)=x1(i) -  (i-1)*(x1(N1)-x1(1))/(N1-1);
0050 end
0051 
0052 % second AR(2) - valley
0053 sigma_e = 0.5;  % process noise
0054 sigma_n = 0.5;  % observation noise
0055 a = [2,-1]';
0056 
0057 x2=zeros(N2,1); % true profile
0058 y2=zeros(N2,1); % observed sample points
0059 for n=2+1:N2
0060     x2(n) = sum(a.*x2(n-1:-1:n-2))+randn(1)*sigma_e;
0061     % discontinuity at position 50
0062     if n==50
0063         x2(n)=x2(n)+5.0;
0064     end
0065     y2(n) = x2(n)+randn(1)*sigma_n;
0066 end
0067 % enforce slope 0
0068 for i=1:N2
0069     y2(i)=y2(i) -  (i-1)*(x2(N2)-x2(1))/(N2-1);
0070     x2(i)=x2(i) -  (i-1)*(x2(N2)-x2(1))/(N2-1);
0071 end
0072 % complete true profile
0073 x = [x1;x2+x1(N1)];
0074 % complete set of observed sample points
0075 y = [y1;y2+y1(N1)];
0076 
0077 % plot the profile
0078 figure('name','Fig. 16.3: Profile Reconstruction','color','w',...
0079     'Position',[0.1*ss(1),0.1*ss(2),0.8*ss(1),0.8*ss(2)]);
0080 
0081 ymin = min(y)-20; ymax = max(y)+20;
0082 subplot(4,3,1); hold on;
0083 plot(1:N,x,'-k','LineWidth',2)
0084 ylim([ymin,ymax])
0085 title('true profile')
0086 
0087 %% sample
0088 disp(' ------ sample profile -----')
0089 Mset = [6,20,48,60,100,144,164,182]';
0090 M = length(Mset);
0091 select = zeros(8,1);
0092 xs = zeros(8,1); ys = zeros(8,1);
0093 for m = 1:M
0094     select(m) = Mset(m);
0095     xs(m) = x(Mset(m));
0096     ys(m) = y(Mset(m));
0097 end
0098 
0099 subplot(4,3,2);hold on
0100 plot(1:N,x,'-k','LineWidth',2)
0101 plot(Mset,ys,'.b','MarkerSize',15)
0102 ylim([ymin,ymax])
0103 title('true profile with sampled points')
0104 
0105 subplot(4,3,3); hold on
0106 plot(Mset,ys,'.b','MarkerSize',15)
0107 ylim([ymin,ymax])
0108 title('sampled points')
0109 
0110 %% interpolation polynomial 0
0111 
0112 disp(' ------ interpolation polynomial 0 -----')
0113 yest = mean(ys);
0114 
0115 subplot(4,3,4); hold on
0116 plot(1:N,yest*ones(N,1),'-r','LineWidth',2)
0117 plot(Mset,ys,'.b','MarkerSize',15)
0118 plot(1:N,x,'--k','LineWidth',1)
0119 ylim([ymin,ymax])
0120 title('best horizontal line')
0121 
0122 %% interpolation polynomial of order 3
0123 
0124 disp(' ------ interpolation order-3-polynomial -----')
0125 order = 3;
0126 % coefficients of sampled coordinates
0127 Xmatrix = [ones(M,1),Mset,Mset.^2,Mset.^3];
0128 % estimate params by regression
0129 estpar  = regress(ys,Xmatrix(:,1:order+1));
0130 
0131 % reconstruction
0132 xc = (1:N)';
0133 Xpred=[ones(N,1),xc,xc.^2,xc.^3];
0134 xpred = Xpred(:,1:order+1)*estpar;
0135 
0136 subplot(4,3,5); hold on
0137 plot(1:N,xpred,'-r','LineWidth',2)   % reconstructed polynomial
0138 plot(Mset,ys,'.b','MarkerSize',15)   % sample points
0139 plot(1:N,x,'--k','LineWidth',1)      % true profile
0140 ylim([ymin,ymax])
0141 title(strcat('best order-',num2str(order),'-polynomial'))
0142 
0143 %% interpolation polynomial of order 6
0144 
0145 disp(' ------ interpolation order-6-polynomial -----')
0146 order = 6;
0147 % coefficients of sampled coordinates
0148 Xmatrix = [ones(M,1),Mset,Mset.^2,Mset.^3,Mset.^4,Mset.^5,Mset.^6];
0149 % estimate params by regression
0150 estpar  = regress(ys,Xmatrix(:,1:order+1));
0151 
0152 % reconstruction
0153 xc=(1:N)';
0154 Xpred=[ones(N,1),xc,xc.^2,xc.^3,xc.^4,xc.^5,xc.^6];
0155 xpred = Xpred(:,1:order+1)*estpar;
0156 
0157 subplot(4,3,6); hold on;
0158 plot(1:N,xpred,'-r','LineWidth',2)
0159 plot(Mset,ys,'.b','MarkerSize',15)
0160 plot(1:N,x,'--k','LineWidth',1)
0161 ylim([ymin,ymax])
0162 title(strcat('best order-',num2str(order),'-polynomial'))
0163 
0164 %% interpolation with very low smoothness
0165 
0166 disp(' ------ interpolation very high smoothness -----')
0167 sigma_e = 0.02;
0168 xest = estimate_profile_robust...
0169     (N,select,ys,sigma_e,1,1,0,[0,0,0,0],0,0);
0170 
0171 subplot(4,3,7);hold on;
0172 plot(1:N,xest,'-r','LineWidth',2)
0173 plot(Mset,ys,'.b','MarkerSize',15)
0174 plot(1:N,x,'--k','LineWidth',1)
0175 ylim([ymin,ymax])
0176 title(['smooth reconstruction, $\sigma_e =',num2str(sigma_e),'$'])
0177 
0178 %% interpolation with medium smoothness
0179 
0180 disp(' ------ interpolation medium smoothness  -----')
0181 sigma_e = 0.1;
0182 xest = estimate_profile_robust...
0183     (N,select,ys,sigma_e,1,1,0,[0,0,0,0],0,0);
0184 
0185 subplot(4,3,8);hold on;
0186 plot(1:N,xest,'-r','LineWidth',2)
0187 plot(Mset,ys,'.b','MarkerSize',15)
0188 plot(1:N,x,'--k','LineWidth',1)
0189 ylim([ymin,ymax])
0190 title(['smooth reconstruction, $\sigma_e =',num2str(sigma_e),'$'])
0191 
0192 %% interpolation with very high smoothness
0193 
0194 disp(' ------ interpolation high smoothness -----')
0195 sigma_e = 0.5;
0196 xest = estimate_profile_robust...
0197     (N,select,ys,sigma_e,1,1,0,[0,0,0,0],0,0);
0198 
0199 subplot(4,3,9);hold on;
0200 plot(1:N,xest,'-r','LineWidth',2)
0201 plot(Mset,ys,'.b','MarkerSize',15)
0202 plot(1:N,x,'--k','LineWidth',1)
0203 ylim([ymin,ymax])
0204 title(['smooth reconstruction, $\sigma_e =',num2str(sigma_e),'$'])
0205 
0206 
0207 %% interpolation medium flatness
0208 
0209 disp(' ------ interpolation high flatness -----')
0210 sigma_e = 0.2;
0211 xest = estimate_profile_robust_flat...
0212     (N,select,ys,sigma_e,1,1,0,[0,0,0,0],0);
0213 
0214 subplot(4,3,10);hold on;
0215 plot(1:N,xest,'-r','LineWidth',2)
0216 plot(Mset,ys,'.b','MarkerSize',15)
0217 plot(1:N,x,'--k','LineWidth',1)
0218 ylim([ymin,ymax])
0219 title(['flat reconstruction, $\sigma_e =',num2str(sigma_e),'$'])
0220 
0221 
0222 %% interpolation with high flatness
0223 
0224 disp(' ------ interpolation low flatness -----')
0225 sigma_e = 1.0;
0226 xest = estimate_profile_robust_flat...
0227     (N,select,ys,sigma_e,1,1,0,[0,0,0,0],0);
0228 
0229 subplot(4,3,11);hold on;
0230 plot(1:N,xest,'-r','LineWidth',2)
0231 plot(Mset,ys,'.b','MarkerSize',15)
0232 plot(1:N,x,'--k','LineWidth',1)
0233 ylim([ymin,ymax])
0234 title(['flat reconstruction, $\sigma_e =',num2str(sigma_e),'$'])
0235 
0236 %% interpolation mixed
0237 
0238 disp(' ------ interpolation left: flat / right: smooth -----')
0239 sigma_e1 = 0.2;
0240 sigma_e2 = 0.5;
0241 fs = [ones(N1,1);ones(N-N1,1)*2];
0242 xest = estimate_profile_robust_flat_smooth...
0243     (N,fs,select,ys,sigma_e1,sigma_e2,1,1,0,[0,0,0,0],0);
0244 
0245 subplot(4,3,12);hold on;
0246 plot(1:N,xest,'-r','LineWidth',2)
0247 plot(Mset,ys,'.b','MarkerSize',15)
0248 plot(1:N,x,'--k','LineWidth',1)
0249 ylim([ymin,ymax])
0250 title(strcat('mixed reconstruction'))
0251 
0252 %% three samples
0253 disp(' --------------------------------------------------')
0254 disp(' --------- Fig. 16.4 Three sample profiles --------')
0255 disp(' --------------------------------------------------')
0256 
0257 N = 200;
0258 N1 = 70;
0259 N2 = N-N1;
0260 
0261 figure('name','Fig. 16.4: Three sample profiles','color','w',...
0262     'Position',[0.3*ss(1),0.2*ss(2),0.5*ss(1),0.3*ss(2)]);
0263 
0264 for samples = 1:3
0265     
0266     [x1,y1] = generate_observed_ARp(N1,1,1,sigma_e,0.5);
0267     % enforce slope 0
0268     for i=1:N1
0269         y1(i)=y1(i) -  (i-1)*(x1(N1)-x1(1))/(N1-1);
0270         x1(i)=x1(i) -  (i-1)*(x1(N1)-x1(1))/(N1-1);
0271     end
0272     
0273     [x2,y2] = generate_observed_ARp(N2,2,1,sigma_n,0.5);
0274     % enforce slope 0
0275     for i=1:N2
0276         y2(i)=y2(i) -  (i-1)*(x2(N2)-x2(1))/(N2-1);
0277         x2(i)=x2(i) -  (i-1)*(x2(N2)-x2(1))/(N2-1);
0278     end
0279     
0280     x = [x1;x2+x1(N1)];
0281     y = [y1;y2+y1(N1)];
0282 
0283     ymin = -250;
0284     ymax = 250;
0285 
0286     subplot(1,3,samples)
0287     hold on
0288     plot(1:N,x,'-k','LineWidth',2)
0289     ylim([ymin,ymax])
0290     title(['sample profile no: ',num2str(samples)])
0291 end