Chapter 5: measuring multiple-source slope parameter

%%%

%%% multi-source slope

%%%

%%% Create DATA

S=3;

T=5; xmin=0; xmax=60; x=linspace(xmin,xmax,T)';

nAlpha=1551; amin=-.025; amax=.025; alpha=[.01 .01 0];

nBeta=101; bmin=97; bmax=100; beta=98.6+.5*(rand(1,S)-.25);

nSig=451; smin=0; smax=.1; SD=.05;

if length(alpha)==1, alpha=alpha*ones(1,S); end; alpha=alpha(1:S);

if length(beta)==1, beta=beta*ones(1,S); end; beta=beta(1:S);

if length(SD)==1, SD=SD*ones(1,S); end; SD=SD(1:S);

noise=nan(T,S);

tmpnoise=rand([T,1]); tmpnoise=tmpnoise-mean(tmpnoise); tmpnoise=tmpnoise/std(tmpnoise);

for s=1:S, noise(:,s)=SD(s)*jumble(tmpnoise,1); end

D=[ones(T,1)*alpha].*[x*ones(1,S)]+[ones(T,1)*beta]+noise;

figure(1); clf; subplot(2,1,1); hold on;

for s=1:S, plot(x,D(:,s),'ko:','MarkerFaceColor',.7*[1 1 1],'MarkerSize',10,'LineWidth',1.25); end

xlim([-10 70]);

tic

L=zeros(nAlpha,nBeta,nSig,S);

Alist=linspace(amin,amax,nAlpha);

Blist=linspace(bmin,bmax,nBeta);

Slist=linspace(smin,smax,nSig+1); Slist=Slist(2:end);

fD=@(Bnow) Bnow+[x*Alist];

for s=1:S+1-(S==1), if s==S+1, Dnow=mean(D,2); else Dnow=D(:,s); end

for nB=1:nBeta, Bnow=Blist(nB);

for nSig=1:nSig, Snow=Slist(nSig);

L(:,nB,nSig,s)=-(T+1)*log(Snow)-.5*(sum((fD(Bnow)-Dnow).^2,1)./Snow^2); end, end, toc, end

if S>1, L1=nan(nAlpha,S);

for snow=1:S,

dprime=alpha(snow)*[0 65]+beta(snow);

plot(x,D(:,snow),'ko:','MarkerFaceColor',.2*[1 1 1],'MarkerSize',10,'LineWidth',2.25);

plot([0 65],dprime,'-','Color',.2*[1 1 1],'LineWidth',2.1);

L1(:,snow)=logsum(L(:,:,:,snow)-max(max(max(L(:,:,:,snow)))),[2 3],[diff(Blist(1:2)) diff(Slist(1:2))]); end

Lavg=logsum(L(:,:,:,end)-max(max(max(L(:,:,:,end)))),[2 3],[diff(Blist(1:2)) diff(Slist(1:2))]);

Lall=sum(logsum(L(:,:,:,1:S)-max(max(max(max(L(:,:,:,1:S))))),[2 3],[diff(Blist(1:2)) diff(Slist(1:2))]),2);

else Lall=logsum(L(:,:,:,1)-max(max(max(L(:,:,:,1)))),[2 3],[diff(Blist(1:2)) diff(Slist(1:2))]); end

figure(2); clf

%amplitude estimate

subplot(2,1,1); hold on

if S>1,

for snow=1:S, plot(Alist,exp(L1(:,snow)-max(L1(:,snow))),'--','Color',[.2 .4 .7]); end

plot(Alist,exp(Lavg-max(Lavg)),':','Color',.7*[1 1 1],'LineWidth',2);

LPVamean=logpeakval([Lavg-max(Lavg) Alist']);

plot(alpha,0,'ko','MarkerFaceColor',.7*[1 1 1]); end

plot(Alist,exp(Lall-max(Lall)),'--','Color',[.2 .4 .7],'LineWidth',2);

LPVa=logpeakval([Lall-max(Lall) Alist']);

plot(alpha,0,'ko','MarkerFaceColor',[.2 .4 .7])

plot([0 0],[0 1],'k--');

axis([-.025 .025 0 1.01])

figure(1); subplot(2,1,2); hold on

dclean=zeros(T,1); Bdelt=mean(D)-mean(x)*LPVa-98.6;

for s=1:S, dclean=dclean+[D(:,s)-Bdelt(s)]; end

plot(x,dclean/S,'ko','MarkerFaceColor',[.2 .3 .6],'MarkerSize',10,'LineWidth',2.25);

plot([-5 0],98.6*[1 1],'-','Color',.5*[1 1 1],'LineWidth',1.2);

dprime=[0 65]*LPVa+98.6; plot([0 65],dprime,'-','Color',.5*[1 1 1],'LineWidth',1.2);

axis([-10 70 98.3 99.5])