% The script file pwsemc04 calculates examples of once and twice
% differentiable processes with dominant frequency, envelopes
% and derivatives.

dt=0.01; eta=1/5; n=300; etadt=eta*dt; gamma=exp(-etadt); 
FI=gamma*[1]; Tp=1; omega=2*pi/Tp;
sigma0=1; 
P=[sigma0*sigma0];
p11=sqrt(P(1,1));
p=[p11];
QQ=P-FI*P*FI';
q11=sqrt(QQ(1,1)); 
Q=[q11]; 
ap=p*[randn]; MA=[ap'];dp=p*[randn];MD=[dp'];

for r=1:n-1
   ab=FI*ap+Q*[randn]; db=FI*dp+Q*[randn];
   MA=[MA;ab']; MD=[MD;dp'];
   ap=ab; dp=db;
end
i=sqrt(-1);
MZ=MA+i*MD;
R=0:n-1;
ME=exp(-i*omega*dt*R);ME=ME';
MZE=MZ(:,1).*ME;
X=real(MZE); Y=imag(MZE); MAB=abs(MZE);
T=0:dt:(n-1)*dt;
plot(T,X,'k',T,Y,'k',T,MAB,'r',T,-MAB,'r')
title('Not differentiable function with dominant frequency')
pause
close




dt=0.01; eta=1/5; n=300; etadt=eta*dt; gamma=exp(-etadt); 
FI=gamma*[1,0;etadt,1]; Tp=1; omega=2*pi/Tp;
sigma1=1; sigma0=sigma1*sqrt(2);
P=[sigma0*sigma0,sigma1*sigma1;sigma1*sigma1,sigma1*sigma1];
p11=sqrt(P(1,1)); p12=0; p21=P(1,2)/p11; p22=sqrt(P(2,2)-p21*p21);
p=[p11,p12;p21,p22];
QQ=P-FI*P*FI';
q11=sqrt(QQ(1,1)); q21=QQ(1,2)/q11; q12=0; q22=sqrt(QQ(2,2)-q21*q21);
Q=[q11,q12;q21,q22]; 
ap=p*[randn;randn]; MA=[ap'];dp=p*[randn;randn];MD=[dp'];
   
for r=1:n-1
   ab=FI*ap+Q*[randn;randn]; db=FI*dp+Q*[randn;randn];
   MA=[MA;ab']; MD=[MD;dp'];
   ap=ab; dp=db;
end
i=sqrt(-1);
MZ=MA+i*MD;
R=0:n-1;
ME=exp(-i*omega*dt*R);ME=ME';
MZE=MZ(:,2).*ME;
X=real(MZE); Y=imag(MZE); MAB=abs(MZE);

T=0:dt:(n-1)*dt;
plot(T,X,'k',T,Y,'k',T,MAB,'r',T,-MAB,'r')
title('Once differentiable function with dominant frequency')
pause
close

DX=(X(2:n,1)-X(1:n-1,1))/dt; DY=(Y(2:n,1)-Y(1:n-1,1))/dt;
DMAB=sqrt(DX.*DX+DY.*DY);
TS=0:dt:(n-2)*dt;TS=TS';
plot(TS,DX,'k',TS,DY,'k',TS,DMAB,'r',TS,-DMAB,'r')
title('Approximate first derivative by finite differeces')
pause
close
plot(TS,DMAB,'r')
title('The envelope of the first derivative')
pause
close


D2X=(X(3:n,1)-2*X(2:n-1,1)+X(1:n-2,1))/dt/dt;
D2Y=(Y(3:n,1)-2*Y(2:n-1,1)+Y(1:n-2,1))/dt/dt;
D2MAB=sqrt(D2X.*D2X+D2Y.*D2Y);
TS2=TS(2:(n-1));
plot(TS2,D2X,'k',TS2,D2Y,'k',TS2,D2MAB,'r',TS2,-D2MAB,'r')
title('Approximate second derivative by finite differeces')
pause
close
plot(TS2,D2MAB,'r.')
title('The envelope of the second derivative')
pause
close



dt=0.01; eta=1/5; n=300; etadt=eta*dt; gamma=exp(-etadt); 
Tp=1; omega=2*pi/Tp;
FI=gamma*[1,0,0;etadt,1,0;etadt*etadt/2,etadt,1]; 
sigma2=1; 
P=sigma2*sigma2*[8/3,4/3,2/3;4/3,4/3,1;2/3,1,1];
p11=sqrt(P(1,1));p12=0;p13=0;p21=P(1,2)/p11;p31=P(1,3)/p11; 
p22=sqrt(P(2,2)-p21*p21);p23=0; p32=(P(2,3)-p31*p21)/p22;
p33=sqrt(P(3,3)-p31*p31-p32*p32);
p=[p11,p12,p13;p21,p22,p23;p31,p32,p33];
QQ=P-FI*P*FI';
q11=sqrt(QQ(1,1)); q21=QQ(1,2)/q11; q31=QQ(1,3)/q11; q12=0; q13=0;
q22=sqrt(QQ(2,2)-q21*q21);q23=0; q32=(QQ(2,3)-q31*q21)/q22;
q33=sqrt(QQ(3,3)-q31*q31-q32*q32);

Q=[q11,q12,q13;q21,q22,q23;q31,q32,q33]; 
ap=p*[randn;randn;randn]; MA=[ap'];
dp=p*[randn;randn;randn]; MD=[dp'];
   
for r=1:n-1
   ab=FI*ap+Q*[randn;randn;randn]; MA=[MA;ab'];
   db=FI*dp+Q*[randn;randn;randn]; MD=[MD;dp'];
   ap=ab; dp=db;
end      
      
i=sqrt(-1);
MZ=MA+i*MD;
R=0:n-1;
ME=exp(-i*omega*dt*R);ME=ME';
MZE=MZ(:,3).*ME;
X=real(MZE); Y=imag(MZE); MAB=abs(MZE);

T=0:dt:(n-1)*dt;
plot(T,X,'k',T,Y,'k',T,MAB,'r',T,-MAB,'r')
title('Twice differentiable function with dominant frequency')
pause
close

DX=(X(2:n,1)-X(1:n-1,1))/dt; DY=(Y(2:n,1)-Y(1:n-1,1))/dt;
DMAB=sqrt(DX.*DX+DY.*DY);
TS=0:dt:(n-2)*dt;TS=TS';
plot(TS,DX,'k',TS,DY,'k',TS,DMAB,'r',TS,-DMAB,'r')
title('Approximate first derivative by finite differeces')
pause
close
plot(TS,DMAB,'r')
title('The envelope of the first derivative')
pause
close

   
D2X=(X(3:n,1)-2*X(2:n-1,1)+X(1:n-2,1))/dt/dt; 
D2Y=(Y(3:n,1)-2*Y(2:n-1,1)+Y(1:n-2,1))/dt/dt;
D2MAB=sqrt(D2X.*D2X+D2Y.*D2Y);
TS2=TS(2:(n-1));
plot(TS2,D2X,'k',TS2,D2Y,'k',TS2,D2MAB,'r',TS2,-D2MAB,'r')
title('Approximate second derivative by finite differeces')
pause
close
plot(TS2,D2MAB,'r')
title('The envelope of the second derivative')
pause
close
