Signal- und Messdatenverarbeitung

Testcodes für Kapitel 7: Kovarianzfunktionen

systematische Zusammenstellung von kombinierbaren Programmabschnitten zur Signal- und Messdatenverarbeitung

Schätzer für verschiedene Signale

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
#python -m pip install --upgrade numpy
#python -m pip install --upgrade matplotlib
#pip install --upgrade numpy
#pip install --upgrade matplotlib
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
%Matlab: erfordert die Statistics and Machine Learning Toolbox
%Octave: pkg load statistics
if exist('OCTAVE_VERSION','builtin')~=0
pkg load statistics
end
periodisches Rauschen mit dem Erwartungswert μ=3 und der Varianz σ2=4
N=1000
x=normal(3,2,N)
dt=1.0
t=arange(0,N)*dt
plot(t,x,'.')
show()
N=1000;
x=normrnd(3,2,1,N);
dt=1.0;
t=(0:N-1)*dt;
plot(t,x,'.')
Korrelationsfunktion
R=zeros(N)
tau=arange(0,N)*dt
for k in range(0,N):
  sum1=0
  for i in range(0,N):
    sum1+=x[i]*x[(i+k)%N]
  
  R[k]=sum1/float(N)

plot(tau,R,'.')
show()
R=zeros(1,N);
tau=(0:N-1)*dt;
for k=0:N-1
sum1=0;
for i=1:N
sum1=sum1+x(i)*x(mod(i+k-1,N)+1);
end
R(k+1)=sum1/N;
end
plot(tau,R,'.')
Kovarianzfunktion
C=zeros(N)
tau=arange(0,N)*dt
mx=mean(x)
for k in range(0,N):
  sum1=0
  for i in range(0,N):
    sum1+=(x[i]-mx)*(x[(i+k)%N]-mx)
  
  C[k]=sum1/float(N)

plot(tau,C,'.')
show()
C=zeros(1,N);
tau=(0:N-1)*dt;
for k=0:N-1
sum1=0;
mx=mean(x);
for i=1:N
sum1=sum1+(x(i)-mx)*(x(mod(i+k-1,N)+1)-mx);
end
C(k+1)=sum1/N;
end
plot(tau,C,'.')
Korrelationskoeffizientenfunktion
rho=zeros(N)
tau=arange(0,N)*dt
mx=mean(x)
for k in range(0,N):
  sum1=0
  for i in range(0,N):
    sum1+=(x[i]-mx)*(x[(i+k)%N]-mx)
  
  rho[k]=sum1/float(N)/var(x)

plot(tau,rho,'.')
show()
rho=zeros(1,N);
tau=(0:N-1)*dt;
for k=0:N-1
sum1=0;
mx=mean(x);
for i=1:N
sum1=sum1+(x(i)-mx)*(x(mod(i+k-1,N)+1)-mx);
end
rho(k+1)=sum1/N/var(x,1);
end
plot(tau,rho,'.')
periodisches Signal (Beispiel harmonisches Signal, mit Rauschen)
N=100
x=sin(2*pi*arange(0,N)/20.0)+normal(0,0.1,N)
dt=1.0
t=arange(0,N)*dt
plot(t,x,'.')
show()
N=100;
x=sin(2*pi*(0:N-1)/20)+normrnd(0,0.1,1,N);
dt=1.0;
t=(0:N-1)*dt;
plot(t,x,'.')
Korrelationsfunktion
R=zeros(N)
tau=arange(0,N)*dt
for k in range(0,N):
  sum1=0
  for i in range(0,N):
    sum1+=x[i]*x[(i+k)%N]
  
  R[k]=sum1/float(N)

plot(tau,R,'.')
show()
R=zeros(1,N);
tau=(0:N-1)*dt;
for k=0:N-1
sum1=0;
for i=1:N
sum1=sum1+x(i)*x(mod(i+k-1,N)+1);
end
R(k+1)=sum1/N;
end
plot(tau,R,'.')
Kovarianzfunktion
C=zeros(N)
tau=arange(0,N)*dt
mx=mean(x)
for k in range(0,N):
  sum1=0
  for i in range(0,N):
    sum1+=(x[i]-mx)*(x[(i+k)%N]-mx)
  
  C[k]=sum1/float(N)

plot(tau,C,'.')
show()
C=zeros(1,N);
tau=(0:N-1)*dt;
for k=0:N-1
sum1=0;
mx=mean(x);
for i=1:N
sum1=sum1+(x(i)-mx)*(x(mod(i+k-1,N)+1)-mx);
end
C(k+1)=sum1/N;
end
plot(tau,C,'.')
Korrelationskoeffizientenfunktion
rho=zeros(N)
tau=arange(0,N)*dt
mx=mean(x)
for k in range(0,N):
  sum1=0
  for i in range(0,N):
    sum1+=(x[i]-mx)*(x[(i+k)%N]-mx)
  
  rho[k]=sum1/float(N)/var(x)

plot(tau,rho,'.')
show()
rho=zeros(1,N);
tau=(0:N-1)*dt;
for k=0:N-1
sum1=0;
mx=mean(x);
for i=1:N
sum1=sum1+(x(i)-mx)*(x(mod(i+k-1,N)+1)-mx);
end
rho(k+1)=sum1/N/var(x,1);
end
plot(tau,rho,'.')
zeitbegrenztes Energiesignal (Beispiel Rechteckimpuls)
N=100
x=zeros(N)
for i in range(10,90):
  x[i]=1

dt=1.0
t=arange(0,N)*dt
plot(t,x,'.')
show()
N=100;
x=zeros(1,N);
for i=10:90
x(i)=1;
end
dt=1.0;
t=(0:N-1)*dt;
plot(t,x,'.')
Korrelationsfunktion
RE=zeros(2*N-1)
tau=zeros(2*N-1)
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0
  for i in range(0,N-abs(k)):
    sum1+=x[i]*x[i+abs(k)]
  
  RE[k]=sum1*dt

plot(tau,RE,'.')
show()
RE=zeros(2*N-1,1);
tau=zeros(2*N-1,1);
for k=-(N-1):N-1
tau(k+N)=k*dt;
sum1=0;
for i=1:N-abs(k)
sum1=sum1+x(i)*x(i+abs(k));
end
RE(k+N)=sum1*dt;
end
plot(tau,RE,'.')
Kovarianzfunktion (für Energiesignale nicht sinnvoll)
CE=zeros(2*N-1)
tau=zeros(2*N-1)
mx=mean(x)
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0
  for i in range(0,N-abs(k)):
    sum1+=(x[i]-mx)*(x[i+abs(k)]-mx)
  
  CE[k]=sum1*dt

plot(tau,CE,'.')
show()
CE=zeros(2*N-1,1);
tau=zeros(2*N-1,1);
mx=mean(x);
for k=-(N-1):N-1
tau(k+N)=k*dt;
sum1=0;
for i=1:N-abs(k)
sum1=sum1+(x(i)-mx)*(x(i+abs(k))-mx);
end
CE(k+N)=sum1*dt;
end
plot(tau,CE,'.')
stationäres Leistungssignal (Beispiel AR1-Prozess)
N=1000
phi=0.95
mx=1
vx=1
theta=sqrt((1-phi**2)*vx)
x=zeros(N)
x[0]=normal(mx,sqrt(vx))
for i in range(1,N):
  x[i]=phi*(x[i-1]-mx)+normal(mx,theta)

dt=1.0
t=arange(0,N)*dt
plot(t,x,'.')
show()
N=1000;
phi=0.95;
mx=1;
vx=1;
theta=sqrt((1-phi^2)*vx);
x=zeros(1,N);
x(1)=normrnd(mx,sqrt(vx));
for i=2:N
x(i)=phi*(x(i-1)-mx)+normrnd(mx,theta);
end
dt=1.0;
t=(0:N-1)*dt;
plot(t,x,'.')
Korrelationsfunktion (erwartungstreu)
R=zeros(2*N-1)
tau=zeros(2*N-1)
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0
  for i in range(0,N-abs(k)):
    sum1+=x[i]*x[i+abs(k)]
  
  R[k]=sum1/(N-abs(k))

plot(tau,R,'.')
show()
R=zeros(2*N-1,1);
tau=zeros(2*N-1,1);
for k=-(N-1):N-1
tau(k+N)=k*dt;
sum1=0;
for i=1:N-abs(k)
sum1=sum1+x(i)*x(i+abs(k));
end
R(k+N)=sum1/(N-abs(k));
end
plot(tau,R,'.')
Kovarianzfunktion (asymptotisch erwartungstreu)
C=zeros(2*N-1)
tau=zeros(2*N-1)
mx=mean(x)
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0
  for i in range(0,N-abs(k)):
    sum1+=(x[i]-mx)*(x[i+abs(k)]-mx)
  
  C[k]=sum1/(N-abs(k))

plot(tau,C,'.')
show()
C=zeros(2*N-1,1);
tau=zeros(2*N-1,1);
mx=mean(x);
for k=-(N-1):N-1
tau(k+N)=k*dt;
sum1=0;
for i=1:N-abs(k)
sum1=sum1+(x(i)-mx)*(x(i+abs(k))-mx);
end
C(k+N)=sum1/(N-abs(k));
end
plot(tau,C,'.')
Korrelationskoeffizientenfunktion (asymptotisch erwartungstreu)
rho=zeros(2*N-1)
tau=zeros(2*N-1)
mx=mean(x)
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0
  for i in range(0,N-abs(k)):
    sum1+=(x[i]-mx)*(x[i+abs(k)]-mx)
  
  rho[k]=sum1/(N-abs(k))/var(x)

plot(tau,rho,'.')
show()
rho=zeros(2*N-1,1);
tau=zeros(2*N-1,1);
mx=mean(x);
for k=-(N-1):N-1
tau(k+N)=k*dt;
sum1=0;
for i=1:N-abs(k)
sum1=sum1+(x(i)-mx)*(x(i+abs(k))-mx);
end
rho(k+N)=sum1/(N-abs(k))/var(x,1);
end
plot(tau,rho,'.')
Korrelationskoeffizientenfunktion (mit lokaler Normierung, asymptotisch erwartungstreu)
rho=zeros(2*N-1)
tau=zeros(2*N-1)
mx=mean(x)
for k in range(-(N-1),N):
  tau[k]=k*dt
  sum1=0
  sum2=0
  sum3=0
  for i in range(0,N-abs(k)):
    sum1+=(x[i]-mx)*(x[i+abs(k)]-mx)
    sum2+=(x[i]-mx)**2
    sum3+=(x[i+abs(k)]-mx)**2
  
  rho[k]=sum1/sqrt(sum2*sum3)

plot(tau,rho,'.')
show()
plot(tau[0:100],rho[0:100],'.',tau[2*N-101:2*N-1],rho[2*N-101:2*N-1],'.')
show()
rho=zeros(2*N-1,1);
tau=zeros(2*N-1,1);
mx=mean(x);
for k=-(N-1):N-1
tau(k+N)=k*dt;
sum1=0;
sum2=0;
sum3=0;
for i=1:N-abs(k)
sum1=sum1+(x(i)-mx)*(x(i+abs(k))-mx);
sum2=sum2+(x(i)-mx)^2;
sum3=sum3+(x(i+abs(k))-mx)^2;
end
rho(k+N)=sum1/sqrt(sum2*sum3);
end
plot(tau,rho,'.')
waitforbuttonpress
plot(tau(N-100:N+100),rho(N-100:N+100),'.')

erwartungstreue Kovarianzschätzung für korrelierte Werte

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
#python -m pip install --upgrade numpy
#python -m pip install --upgrade matplotlib
#pip install --upgrade numpy
#pip install --upgrade matplotlib
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
%Octave: pkg load statistics
%Matlab: erfordert die Statistics and Machine Learning Toolbox
Generieren einer korrelierten Zeitreihe (Moving-Average-Prozess)
N=50
Q=20
mu=3
sigma=1
e=normal(0,sigma/sqrt(Q),int(N+Q))
x=zeros(N)
x[0]=sum(e[0:Q])+mu
for i in range(1,N):
  x[i]=x[i-1]-e[i-1]+e[i-1+Q]

plot(x,'o')
show()
N=50;
Q=20;
mu=3;
sigma=1;
e=normrnd(0,sigma/sqrt(Q),[1,N+Q]);
x=zeros(1,N);
x(1)=sum(e(1:Q))+mu;
for i=2:N
x(i)=x(i-1)-e(i-1)+e(i-1+Q);
end
plot(x,'o')
Bestimmung der primären Kovarianzfunktion
K=60
mest=mean(x)
tau=zeros(K)
Cprim=zeros(K)
Csim=zeros(K)
for k in range(-(K//2),(K+1)//2):
  tau[k]=k
  Csim[k]=maximum(0,1-abs(tau[k])/float(Q))
  for i in range(0,N-abs(k)):
    Cprim[k]+=(x[i]-mest)*(x[i+abs(k)]-mest)/(N-abs(k))
  

plot(roll(tau,(K+1)//2),roll(Csim,(K+1)//2),tau,Cprim,'o')
show()
K=60;
mest=mean(x);
tau=zeros(1,K);
Cprim=zeros(1,K);
for k=-fix(K/2):K-1-fix(K/2)
tau(mod(k,K)+1)=k;
for i=1:N-abs(k)
Cprim(mod(k,K)+1)=Cprim(mod(k,K)+1)+(x(i)-mest)*(x(i+abs(k))-mest)/(N-abs(k));
end
end
plot(circshift(tau,[0;fix(K/2)]),circshift(max(0,1-abs(tau)/Q),[0;fix(K/2)]),tau,Cprim,'o')
Bestimmung der Abbildungsmatrix und deren Inversen
from numpy.linalg import *
A=zeros([K,K])
for i in range(-(K//2),(K+1)//2):
  for j in range(-(K//2),(K+1)//2):
    A[i,j]=int(i==j)-2*(N-maximum(maximum(abs(i),abs(j)),minimum(N,abs(i-j))))/float(N*(N-abs(i)))+(N-abs(j))/float(N**2)
  

Ainv=inv(A)
A=zeros(K,K);
for i=-fix(K/2):K-1-fix(K/2)
for j=-fix(K/2):K-1-fix(K/2)
A(mod(i,K)+1,mod(j,K)+1)=(i==j)-2*(N-max(max(abs(i),abs(j)),min(N,abs(i-j))))/(N*(N-abs(i)))+(N-abs(j))/N^2;
end
end
Ainv=inv(A);
Korrektur der Kovarianzfunktion
Ccorr=dot(Ainv,Cprim)
plot(roll(tau,(K+1)//2),roll(Csim,(K+1)//2),tau,Cprim,'o',tau,Ccorr,'o')
show()
Ccorr=Ainv*Cprim';
plot(circshift(tau,[0;fix(K/2)]),circshift(max(0,1-abs(tau)/Q),[0;fix(K/2)]),tau,Cprim,'o',tau,Ccorr,'o')
Nachweis der Erwartungstreue der Korrektur
from numpy.linalg import *
M=1000
N=50
Q=20
K=60
mu=3
sigma=1
tau=zeros(K)
Csim=zeros(K)
for k in range(-(K//2),(K+1)//2):
  tau[k]=k
  Csim[k]=maximum(0,1-abs(tau[k])/float(Q))

A=zeros([K,K])
for k in range(-(K//2),(K+1)//2):
  for j in range(-(K//2),(K+1)//2):
    A[k,j]=int(k==j)+(N-abs(j))/float(N**2)-(minimum(minimum(N,N-j),N-k)-maximum(maximum(0,-j),-k)+minimum(minimum(N,N-j),N+k-j)-maximum(maximum(0,-j),k-j))/float(N*(N-abs(k)))
  

Ainv=inv(A)
SCprim=zeros(K)
SCcorr=zeros(K)
x=zeros(N)
for m in range(0,M):
  e=normal(0,sigma/sqrt(Q),int(N+Q))
  x[0]=sum(e[0:Q])+mu
  for i in range(1,N):
    x[i]=x[i-1]-e[i-1]+e[i-1+Q]
  
  mest=mean(x)
  Cprim=zeros(K)
  for k in range(-(K//2),(K+1)//2):
    for i in range(0,N-abs(k)):
      Cprim[k]+=(x[i]-mest)*(x[i+abs(k)]-mest)/(N-abs(k))
    
  
  SCprim+=Cprim
  Ccorr=dot(Ainv,Cprim)
  SCcorr+=Ccorr

plot(roll(tau,(K+1)//2),roll(Csim,(K+1)//2),tau,SCprim/float(M),'o',tau,SCcorr/float(M),'o')
show()
M=1000;
N=50;
Q=20;
K=60;
mu=3;
sigma=1;
tau=zeros(1,K);
for k=-fix(K/2):K-1-fix(K/2)
tau(mod(k,K)+1)=k;
end
A=zeros(K,K);
for k=-fix(K/2):K-1-fix(K/2)
for j=-fix(K/2):K-1-fix(K/2)
A(mod(k,K)+1,mod(j,K)+1)=(k==j)+(N-abs(j))/N^2-(min(min(N,N-j),N-k)-max(max(0,-j),-k)+min(min(N,N-j),N+k-j)-max(max(0,-j),k-j))/(N*(N-abs(k)));
end
end
Ainv=inv(A);
SCprim=zeros(1,K);
SCcorr=zeros(1,K);
x=zeros(1,N);
for m=1:M
e=normrnd(0,sigma/sqrt(Q),[1,N+Q]);
x(1)=sum(e(1:Q))+mu;
for i=2:N
x(i)=x(i-1)-e(i-1)+e(i-1+Q);
end
mest=mean(x);
Cprim=zeros(1,K);
for k=-fix(K/2):K-1-fix(K/2)
for i=1:length(x)-abs(k)
Cprim(mod(k,K)+1)=Cprim(mod(k,K)+1)+(x(i)-mest)*(x(i+abs(k))-mest)/(N-abs(k));
end
end
SCprim=SCprim+Cprim;
Ccorr=Ainv*Cprim';
SCcorr=SCcorr+Ccorr;
end
plot(circshift(tau,[0;fix(K/2)]),circshift(max(0,1-abs(tau)/Q),[0;fix(K/2)]),tau,SCprim/M,'o',tau,SCcorr/M,'o')
Vergleich des mittleren quadratischen Fehlers
from numpy.linalg import *
M=1000
N=50
Q=20
K=60
mu=3
sigma=1
tau=zeros(K)
Csim=zeros(K)
for k in range(-(K//2),(K+1)//2):
  tau[k]=k
  Csim[k]=maximum(0,1-abs(tau[k])/float(Q))

A=zeros([K,K])
for k in range(-(K//2),(K+1)//2):
  for j in range(-(K//2),(K+1)//2):
    A[k,j]=int(k==j)+(N-abs(j))/float(N**2)-(minimum(minimum(N,N-j),N-k)-maximum(maximum(0,-j),-k)+minimum(minimum(N,N-j),N+k-j)-maximum(maximum(0,-j),k-j))/float(N*(N-abs(k)))
  

Ainv=inv(A)
SCprim=zeros(K)
SCcorr=zeros(K)
x=zeros(N)
for m in range(0,M):
  e=normal(0,sigma/sqrt(Q),int(N+Q))
  x[0]=sum(e[0:Q])+mu
  for i in range(1,N):
    x[i]=x[i-1]-e[i-1]+e[i-1+Q]
  
  mest=mean(x)
  Cprim=zeros(K)
  for k in range(-(K//2),(K+1)//2):
    for i in range(0,N-abs(k)):
      Cprim[k]+=(x[i]-mest)*(x[i+abs(k)]-mest)/(N-abs(k))
    
  
  SCprim+=(Cprim-Csim)**2
  Ccorr=dot(Ainv,Cprim)
  SCcorr+=(Ccorr-Csim)**2

plot(tau,SCprim/float(M),'o',tau,SCcorr/float(M),'x')
show()
M=1000;
N=50;
Q=20;
K=60;
mu=3;
sigma=1;
tau=zeros(1,K);
for k=-fix(K/2):K-1-fix(K/2)
tau(mod(k,K)+1)=k;
end
A=zeros(K,K);
for k=-fix(K/2):K-1-fix(K/2)
for j=-fix(K/2):K-1-fix(K/2)
A(mod(k,K)+1,mod(j,K)+1)=(k==j)+(N-abs(j))/N^2-(min(min(N,N-j),N-k)-max(max(0,-j),-k)+min(min(N,N-j),N+k-j)-max(max(0,-j),k-j))/(N*(N-abs(k)));
end
end
Ainv=inv(A);
SCprim=zeros(1,K);
SCcorr=zeros(1,K);
x=zeros(1,N);
for m=1:M
e=normrnd(0,sigma/sqrt(Q),[1,N+Q]);
x(1)=sum(e(1:Q))+mu;
for i=2:N
x(i)=x(i-1)-e(i-1)+e(i-1+Q);
end
mest=mean(x);
Cprim=zeros(1,K);
for k=-fix(K/2):K-1-fix(K/2)
for i=1:length(x)-abs(k)
Cprim(mod(k,K)+1)=Cprim(mod(k,K)+1)+(x(i)-mest)*(x(i+abs(k))-mest)/(N-abs(k));
end
end
SCprim=SCprim+(Cprim-max(0,1-abs(tau)/Q)).^2;
Ccorr=(Ainv*Cprim')';
SCcorr=SCcorr+(Ccorr-max(0,1-abs(tau)/Q)).^2;
end
plot(tau,SCprim/M,'o',tau,SCcorr/M,'x')