Signal- und Messdatenverarbeitung

Lineare stochastische Prozesse

ARp-Prozess

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
%Octave: pkg load statistics
%Matlab: erfordert die Statistics and Machine Learning Toolbox
Vorgabe einer Korrelationskoeffizientenfunktion (einseitig)
dt=1.0
K=20
rho=zeros(K)
for k in range(0,K):
  rho[k]=exp(-(k/5.0)**2-k/100.0) # rho[0]=1

plot(arange(0,K)*dt,rho,'o')
show()
dt=1;
K=20;
rho=zeros(1,K);
for k=1:K
rho(k)=exp(-((k-1)/5)^2-(k-1)/100); % rho(1)=1
end
plot((0:K-1)*dt,rho,'o')
Für eine Prozessidentifikation aus einer Zeitreihe von Werten sind die Korrelationskoeffizientenfunktion und die Varianz aus dem Datensatz zu schätzen (für reelle Signale ohne Gewichtung bzw. mit Gewichtung).
Entwurf der Koeffizienten eines ARp-Prozesses aus der vorgegebenen Korrelationskoeffizientenfunktion sowie der Varianz
from numpy.linalg import *
rho/=float(rho[0]) # rho[0]=1
v=1.0 # Varianz
p=len(rho)-1
RHO=zeros([p,p])
for i in range(0,p):
  for j in range(0,p):
    RHO[i,j]=rho[abs(i-j)];
  

eq=zeros(p)
for i in range(0,p):
  eq[i]=rho[i+1];

phi=real(linalg.solve(RHO,eq))
theta0=real(sqrt(v*(1-sum(phi*rho[1:p+1]))))
plot(phi,'o')
show()
rho=rho/rho(1); % rho(1)=1
v=1; % Varianz
p=length(rho)-1;
RHO=zeros(p,p);
for i=1:p
for j=1:p
RHO(i,j)=rho(abs(i-j)+1);
end
end
eq=zeros(1,p);
for i=1:p
eq(i)=rho(i+1);
end
phi=linsolve(RHO,eq')';
theta0=sqrt(v*(1-sum(phi.*rho(2:p+1))));
plot(phi,'o')
Einzelrealisierung (inklusive der Initialisierung)
from numpy.linalg import *
p=len(phi)
PHI=zeros([p,p])
for i in range(0,p):
  for j in range(0,p):
    PHI[i,j]=0
    if (i+j+1>=0) and (i+j+1<p):
      PHI[i,j]+=phi[i+j+1]
    
    if (i-j-1>=0) and (i-j-1<p):
      PHI[i,j]+=phi[i-j-1]
    
    if i==j:
      PHI[i,j]-=1
    
  

eq=zeros(p)
for i in range(0,p):
  eq[i]=-phi[i]

rho=append([1],linalg.solve(PHI,eq))
RHO=zeros([p,p])
for i in range(0,p):
  for j in range(0,p):
    RHO[i,j]=rho[abs(i-j)];
  

A=zeros([p,p])
for i in range(0,p):
  A[i,i]=1.0

bl=inf
B=(dot(inv(A.T),RHO)-A)/2.0
b=amax(abs(B))
A=A+B
while b<bl:
  bl=b
  B=(dot(inv(A.T),RHO)-A)/2.0
  b=amax(abs(B))
  A=A+B

e=standard_normal(p)
v=theta0**2/float(1-sum(phi*rho[1:p+1]))
u=sqrt(v)*dot(A,e.T)

N=1000
m=3.0 # Erwartungswert
x=zeros(N)
for i in range(0,N):
  u=append(sum(u*phi)+normal(0,theta0),u[0:p-1])
  x[i]=u[0]+m

t=arange(0,N)*dt
plot(t,x,'o')
show()
p=length(phi);
PHI=zeros(p,p);
for i=1:p
for j=1:p
PHI(i,j)=0;
if (i+j>0) && (i+j<=p)
PHI(i,j)=PHI(i,j)+phi(i+j);
end
if (i-j>0) && (i-j<=p)
PHI(i,j)=PHI(i,j)+phi(i-j);
end
if i==j
PHI(i,j)=PHI(i,j)-1;
end
end
end
eq=zeros(1,p);
for i=1:p
eq(i)=-phi(i);
end
rho=[1,linsolve(PHI,eq')'];
RHO=zeros(p,p);
for i=1:p
for j=1:p
RHO(i,j)=rho(abs(i-j)+1);
end
end
A=sqrtm(RHO);
e=randn(1,p);
v=theta0^2/(1-sum(phi.*rho(2:p+1)));
u=(sqrt(v)*A*e')';
N=1000;
m=3 % Erwartungswert
x=zeros(1,N);
for i=1:N
u=[sum(u.*phi)+normrnd(0,theta0),u(1:p-1)];
x(i)=u(1)+m;
end
t=(0:N-1)*dt;
plot(t,x,'o')
Kovarianzfunktion des entworfenen Prozesses
from numpy.linalg import *
p=len(phi)
PHI=zeros([p,p])
for i in range(0,p):
  for j in range(0,p):
    PHI[i,j]=0
    if (i+j+1>=0) and (i+j+1<p):
      PHI[i,j]+=phi[i+j+1]
    
    if (i-j-1>=0) and (i-j-1<p):
      PHI[i,j]+=phi[i-j-1]
    
    if i==j:
      PHI[i,j]-=1
    
  

eq=zeros(p)
for i in range(0,p):
  eq[i]=-phi[i]

rho=append([1],linalg.solve(PHI,eq))
v=theta0**2/float(1-sum(phi*rho[1:p+1]))
K=100
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
C=zeros(K)
for k in range(maximum(-(K//2),-p),minimum((K+1)//2,p+1)):
  C[k]=v*rho[abs(k)]

for k in range(minimum((K+1)//2,p+1),(K+1)//2):
  C[k]=0
  for j in range(0,p):
    C[k]+=phi[j]*C[k-(j+1)]
  

for k in range(maximum(-(K//2),-p-1),-(K//2)-1,-1):
  C[k]=0
  for j in range(0,p):
    C[k]+=phi[j]*C[k+j+1]
  

plot(tau,C,'o')
show()
p=length(phi);
PHI=zeros(p,p);
for i=1:p
for j=1:p
PHI(i,j)=0;
if (i+j>0) && (i+j<=p)
PHI(i,j)=PHI(i,j)+phi(i+j);
end
if (i-j>0) && (i-j<=p)
PHI(i,j)=PHI(i,j)+phi(i-j);
end
if i==j
PHI(i,j)=PHI(i,j)-1;
end
end
end
eq=zeros(1,p);
for i=1:p
eq(i)=-phi(i);
end
rho=[1,linsolve(PHI,eq')'];
v=theta0^2/(1-sum(phi.*rho(2:p+1)));
K=100;
tau=circshift((-fix(K/2):K-1-fix(K/2))*dt,[0;fix((K+1)/2)]);
C=zeros(1,K);
for k=max(-fix(K/2),-p):min(K-1-fix(K/2),p)
C(mod(k,K)+1)=v*rho(abs(k)+1);
end
for k=min(K-fix(K/2),p+1):K-1-fix(K/2)
C(mod(k,K)+1)=0;
for j=1:p
C(mod(k,K)+1)=C(mod(k,K)+1)+phi(j)*C(mod(k-j,K)+1);
end
end
for k=max(-fix(K/2),-p)-1:-1:-fix(K/2)
C(mod(k,K)+1)=0;
for j=1:p
C(mod(k,K)+1)=C(mod(k,K)+1)+phi(j)*C(mod(k+j,K)+1);
end
end
plot(tau,C,'o')
Leistungsdichtespektrums des entworfenen Prozesses
K=1000
p=len(phi)
df=1/float(K*dt)
f=roll(arange(-(K//2),(K+1)//2)*df,(K+1)//2)
S=zeros(K)
for i in range(0,len(f)):
  S[i]=dt*theta0**2/abs(1-sum(phi*exp(-2j*pi*arange(1,p+1)*f[i]*dt)))**2

plot(f,S,'o')
show()
K=1000;
p=length(phi);
df=1/(K*dt);
f=circshift((-fix(K/2):fix((K-1)/2))*df,[0;fix((K+1)/2)]);
S=zeros(1,K);
for i=(1:length(f))
S(i)=dt*theta0^2/abs(1-sum(phi.*exp(-2j*pi*(1:p)*f(i)*dt)))^2;
end
plot(f,S,'o')

MAq-Prozess

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
Vorgabe des Leistungsdichtespektrums (beliebige, paarweise verschiedene Frequenzen)
dt=1.0
J=20
df=0.024
f=arange(1,J+1)*df
S=0.1/(0.02+f**(5.0/3.0))*exp(-10*f)
plot(f,S,'o')
show()
dt=1.0;
J=20;
df=0.024;
f=(1:J)*df;
S=0.1./(0.02+f.^(5.0/3.0)).*exp(-10*f);
plot(f,S,'o')
Entwurf der Koeffizienten eines MAq-Prozesses aus der vorgegebenen Leistungsdichte
COS=zeros([J,J])
for i in range(0,J):
  COS[i,0]=1
  for j in range(1,J):
    COS[i,j]=2*cos(2*pi*f[i]*j*dt)
  

eq=sqrt(S/dt)
q=2*J-2
theta=zeros(q+1)
theta[J-1:q+1]=linalg.solve(COS,eq)
for i in range(0,q+1):
  theta[i]=theta[abs(q-i)]

plot(theta,'o')
show()
COS=zeros(J,J);
for i=1:J
COS(i,1)=1;
for j=2:J
COS(i,j)=2*cos(2*pi*f(i)*(j-1)*dt);
end
end
eq=sqrt(S/dt);
q=2*J-2;
theta=zeros(1,q+1);
theta(J:q+1)=linsolve(COS,eq')';
for i=1:q+1
theta(i)=theta(abs(q-i)+2);
end
plot(theta,'o')
Vorgabe der Kovarianzfunktion
dt=1.0
K=20
tau=arange(0,K)*dt
C=2*(1-arange(0,K)/float(K))*cos(1.5*pi*arange(0,K)/float(K))
plot(tau,C,'o')
show()
dt=1.0;
K=20;
tau=(0:K-1)*dt;
C=2*(1-(0:K-1)/K).*cos(1.5*pi*(0:K-1)/K);
plot(tau,C,'o')
Entwurf der Koeffizienten eines MAq-Prozesses aus der vorgegebenen Kovarianzfunktion
from numpy.linalg import *
q=K-1
theta=sqrt(C[0]/float(q+1))*(q+1-arange(0,q+1))/float(q+1)
for l in range(0,20):
  J=zeros([q+1,q+1])
  e=-C[0:q+1]
  for i in range(0,q+1):
    for j in range(i,q+1):
      e[i]+=theta[j-i]*theta[j]
    for j in range(0,q+1):
      if (j-i>=0) and (j-i<=q):
        J[i,j]+=theta[j-i]
      if (j+i>=0) and (j+i<=q):
        J[i,j]+=theta[j+i]
      
    
  theta+=-dot(inv(J),e)

plot(theta,'o')
show()
q=K-1;
theta=sqrt(C(1)/(q+1))*ones(1,q+1);
Theta=zeros(q+1,q+1);
for l=1:1000
for i=1:q+1
for j=i:q+1
Theta(i,j)=theta(j-i+1);
end
end
theta=(2*theta+linsolve(Theta,C')')/3;
end
plot(theta,'o')
Einzelrealisierung (inklusive der Initialisierung)
from numpy.fft import *
N=1000
m=3 # Erwartungswert
q=len(theta)-1
e=standard_normal(N+q+1)
xp=real(ifft(fft(e)*fft(append(theta,zeros(N)))))
x=xp[0:N]
t=arange(0,N)*dt
plot(t,x,'o')
show()
N=1000;
m=3; % Erwartungswert
q=length(theta)-1;
e=randn(1,N+q+1);
xp=real(ifft(fft(e).*fft([theta,zeros(1,N)])));
x=xp(1:N)+m;
t=(0:N-1)*dt;
plot(t,x,'o')
Kovarianzfunktion des entworfenen Prozesses
K=100
q=len(theta)-1
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
C=zeros(K)
for k in range(-(K//2),(K+1)//2):
  C[k]=0
  for i in range(0,q+1-abs(k)):
    C[k]+=theta[i]*theta[i+abs(k)]
  

plot(tau,C,'o')
show()
K=100;
q=length(theta)-1;
tau=circshift((-fix(K/2):K-1-fix(K/2))*dt,[0;fix((K+1)/2)]);
C=zeros(1,K);
for k=-fix(K/2):K-fix(K/2)
C(mod(k,K)+1)=0;
for i=0:q-abs(k)
C(mod(k,K)+1)=C(mod(k,K)+1)+theta(i+1)*theta(i+1+abs(k));
end
end
plot(tau,C,'o')
Leistungsdichtespektrum des entworfenen Prozesses
K=1000
q=len(theta)-1
df=1/float(K*dt)
f=roll(arange(-(K//2),(K+1)//2)*df,(K+1)//2)
S=zeros(K)
for i in range(0,len(f)):
  S[i]=dt*abs(sum(theta*exp(-2j*pi*arange(0,q+1)*f[i]*dt)))**2

plot(f,S,'o')
show()
K=1000;
q=length(theta)-1;
df=1/(K*dt);
f=circshift((-fix(K/2):fix((K-1)/2))*df,[0;fix((K+1)/2)]);
S=zeros(1,K);
for i=(1:length(f))
S(i)=dt*abs(sum(theta.*exp(-2j*pi*(0:q)*f(i)*dt)))^2;
end
plot(f,S,'o')