Signal- und Messdatenverarbeitung

Testcodes für Kapitel 21: Ungleichmäßige Abtastung

systematische Zusammenstellung von kombinierbaren Programmabschnitten zur Signal- und Messdatenverarbeitung

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 numpy.fft 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
Generieren von Wertepaaren (ti,xi) mit den Messzeitpunkten ti und den Messwerten xi (AR1-Prozess als Beispiel, mit Mittelwert verschieden von null, kurzer Datensatz zur Visualisierung)
T=100.0
Ti=10.0
dd=1.0
mxs=3.0
vxs=1.0
t=[]
x=[]
te=0.0
xe=normal(mxs,sqrt(vxs))
while te<T:
  tp=exponential(1.0/dd)
  te+=tp
  if te<T:
    t.append(te)
    phi=exp(-tp/Ti)
    theta=sqrt((1-phi**2)*vxs)
    xe=(xe-mxs)*phi+normal(mxs,theta)
    x.append(xe)
  

t=array(t)
x=array(x)
plot(t,x,'o',t,zeros(len(t)),'.')
show()
T=100;
Ti=10;
dd=1;
mxs=3;
vxs=1;
t=[];
x=[];
te=0;
xe=normrnd(mxs,sqrt(vxs));
while te<T
tp=exprnd(1/dd);
te=te+tp;
if te<T
t(end+1)=te;
phi=exp(-tp/Ti);
theta=sqrt((1-phi^2)*vxs);
xe=(xe-mxs)*phi+normrnd(mxs,theta);
x(end+1)=xe;
end
end
plot(t,x,'o',t,zeros(1,length(t)),'.')
Generieren von Wertepaaren (ti,xi) mit den Messzeitpunkten ti und den Messwerten xi (AR1-Prozess als Beispiel, mit Mittelwert verschieden von null, langer Datensatz zur Weiterverarbeitung)
T=10000.0
Ti=10.0
dd=1.0
mxs=3.0
vxs=1.0
t=[]
x=[]
te=0.0
xe=normal(mxs,sqrt(vxs))
while te<T:
  tp=exponential(1.0/dd)
  te+=tp
  if te<T:
    t.append(te)
    phi=exp(-tp/Ti)
    theta=sqrt((1-phi**2)*vxs)
    xe=(xe-mxs)*phi+normal(mxs,theta)
    x.append(xe)
  

t=array(t)
x=array(x)
plot(t,x,'.')
show()
T=10000;
Ti=10;
dd=1;
mxs=3;
vxs=1;
t=[];
x=[];
te=0;
xe=normrnd(mxs,sqrt(vxs));
while te<T
tp=exprnd(1/dd);
te=te+tp;
if te<T
t(end+1)=te;
phi=exp(-tp/Ti);
theta=sqrt((1-phi^2)*vxs);
xe=(xe-mxs)*phi+normrnd(mxs,theta);
x(end+1)=xe;
end
end
plot(t,x,'.')
Mittelwert x‾=1N⋅∑i=0N−1xi
mxe=mean(x)
print(mxe,mxs)
mxe=mean(x);
[mxe,mxs]
Varianz (ohne Bessel-Korrektur, asymptotisch erwartungstreu) s2=1N⋅∑i=0N−1(xi−x‾)2=(1N⋅∑i=0N−1xi2)−x‾2
vxe=var(x)
print(vxe,vxs)
vxe=var(x,1);
[vxe,vxs]
Autokorrelationsfunktion und Leistungsdichtespektrum über Slotkorrelation (ohne Selbstprodukte durch bk(0)=0) Rk=R(τk)=∑i=0N−1∑j=0N−1bk(tj−ti)xixj∑i=0N−1∑j=0N−1bk(tj−ti)
bk(t)={1falls 0<|t|<Δt/20sonst
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K<2T/Δt
Sj=S(fj)=Δt⋅FFT{Rk}=Δt⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Rk𝐞−2π𝐢jk/K
(imaginäre Einheit 𝐢)
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
N=len(t)
dt=1.0
K=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R1=zeros(K)
R0=zeros(K)
i=1
while i<N:
  j=i-1
  while (j>=0) and (t[j]>t[i]-(K//2+0.5)*dt):
    k=int(round((t[j]-t[i])/float(dt)))
    R1[k]+=x[i]*x[j]
    R0[k]+=1.0
    j-=1
  
  i+=1

for k in range(1,(K+1)//2):
  R1[k]=R1[-k]
  R0[k]=R0[-k]

R=R1/R0
plot(tau,R,'.',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs**2)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'.',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
K=200;
K=min(ceil(2*T/dt)-1,K);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R1=zeros(1,K);
R0=zeros(1,K);
i=2;
while i<=N
j=i-1;
while (j>0) && (t(j)>t(i)-(fix(K/2)+0.5)*dt)
k=round((t(j)-t(i))/dt);
R1(mod(K+k,K)+1)=R1(mod(K+k,K)+1)+x(i)*x(j);
R0(mod(K+k,K)+1)=R0(mod(K+k,K)+1)+1;
j=j-1;
end
i=i+1;
end
for k=2:fix((K+1)/2)
R1(k)=R1(K-k+2);
R0(k)=R0(K-k+2);
end
R=R1./R0;
plot(tau,R,'.',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs^2)
pause
%
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'.',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, ohne Abzug der Selbstprodukte, ohne Normierung) Ej=E(fj)=T2N(N−1)|Xj|2
Xj=∑i=0N−1xi𝐞−2π𝐢fjti
(imaginäre Einheit 𝐢)
fj=jJΔt
j=−⌊J/2⌋…⌊(J−1)/2⌋
J≥⌈2T/Δt⌉
RE,k=RE(τk)=1Δt⋅IFFT{Ej}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Ej𝐞2π𝐢fjτk/J
Rk=R(τk)=RE,kT−|τk|
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K<2T/Δt
Sj=S(fj)=Δt⋅FFT{Rk}=Δt⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Rk𝐞−2π𝐢jk/K
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
N=len(t)
dt=1.0
J=int(ceil(2*T/float(dt)))
fp=arange(0,J//2+1)/float(J*dt)
X=zeros(size(fp))+0j
for i in range(0,N):
  X+=x[i]*exp(-2j*pi*fp*t[i])

E=zeros(J)
for k in range(0,J//2+1):
  E[k]=T**2*(abs(X[k])**2)/float(N**2)

for k in range(-(J//2-1),0):
  E[k]=E[-k]

RE=real(ifft(E))/float(dt)
K=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)
for k in range(-(K//2),(K+1)//2):
  R[k]=RE[k]/float(T-abs(tau[k]))

plot(tau,R,'.',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs**2)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'.',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
J=ceil(2*T/dt);
fp=(0:fix(J/2))/(J*dt);
X=zeros(size(fp))+0j;
for i=1:N
X=X+x(i)*exp(-2j*pi*fp*t(i));
end
E=zeros(1,J);
for k=1:fix(J/2)+1
E(k)=T^2*(abs(X(k)).^2)/(N^2);
end
for k=fix(J/2)+2:J
E(k)=E(J-k+2);
end
RE=real(ifft(E))/dt;
K=200;
K=min(ceil(2*T/dt)-1,K);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros(1,K);
for k=1:fix(K/2)+1
R(k)=RE(k)/(T-abs(tau(k)));
end
for k=fix(K/2)+2:K
R(k)=R(K-k+2);
end
plot(tau,R,'.',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs^2)
pause
%
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'.',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Fourier-Transformation, mit Abzug der Selbstprodukte, ohne Normierung) Ej=E(fj)=T2N(N−1){|Xj|2−∑i=0N−1xi2}
Xj=∑i=0N−1xi𝐞−2π𝐢fjti
(imaginäre Einheit 𝐢)
fj=jJΔt
j=−⌊J/2⌋…⌊(J−1)/2⌋
J≥⌈2T/Δt⌉
RE,k=RE(τk)=1Δt⋅IFFT{Ej}=1JΔt⋅∑j=−⌊J/2⌋⌊(J−1)/2⌋Ej𝐞2π𝐢fjτk/J
Rk=R(τk)=RE,kT−|τk|
τk=kΔt
k=−⌊K/2⌋…⌊(K−1)/2⌋
K<2T/Δt
Sj=S(fj)=Δt⋅FFT{Rk}=Δt⋅∑k=⌊K/2⌋⌊(K−1)/2⌋Rk𝐞−2π𝐢jk/K
fj=jKΔt
j=−⌊K/2⌋…⌊(K−1)/2⌋
N=len(t)
dt=1.0
J=int(ceil(2*T/float(dt)))
fp=arange(0,J//2+1)/float(J*dt)
X=zeros(size(fp))+0j
for i in range(0,N):
  X+=x[i]*exp(-2j*pi*fp*t[i])

E=zeros(J)
for k in range(0,J//2+1):
  E[k]=T**2*(abs(X[k])**2-sum(x**2))/float(N**2-N)

for k in range(-(J//2-1),0):
  E[k]=E[-k]

RE=real(ifft(E))/float(dt)
K=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)
for k in range(-(K//2),(K+1)//2):
  R[k]=RE[k]/float(T-abs(tau[k]))

plot(tau,R,'.',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs**2)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'.',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
J=ceil(2*T/dt);
fp=(0:fix(J/2))/(J*dt);
X=zeros(size(fp))+0j;
for i=1:N
X=X+x(i)*exp(-2j*pi*fp*t(i));
end
E=zeros(1,J);
for k=1:fix(J/2)+1
E(k)=T^2*(abs(X(k)).^2-sum(x.^2))/(N^2-N);
end
for k=fix(J/2)+2:J
E(k)=E(J-k+2);
end
RE=real(ifft(E))/dt;
K=200;
K=min(ceil(2*T/dt)-1,K);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros(1,K);
for k=1:fix(K/2)+1
R(k)=RE(k)/(T-abs(tau(k)));
end
for k=fix(K/2)+2:K
R(k)=R(K-k+2);
end
plot(tau,R,'.',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs^2)
pause
%
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'.',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Lomb-Scargle-Methode, ohne Abzug der Selbstprodukte, ohne Normierung, ohne Korrektur des dynamischen Fehlers)
from numpy.fft import *
import scipy.signal as scisig
N=len(x)
dt=1.0
mxe=mean(x)
J=int(ceil(2*T/float(dt)))
fp=arange(1,J//2+1)/float(J*dt)
X=scisig.lombscargle(t,x-mxe,2*pi*fp)
E=zeros(J)
for k in range(0,J//2+1):
  E[k]=T**2*(N*X[k-1])/float(N**2)

for k in range(-((J-1)//2),0):
  E[k]=E[-k]

RE=real(ifft(E))/float(dt)
K=200
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)
for k in range(-(K//2),(K+1)//2):
  R[k]=RE[k]/float(T-abs(tau[k]))+mxe**2

plot(tau,R,'.',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs**2)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'.',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(x);
dt=1.0;
mxe=mean(x);
J=ceil(2*T/dt);
fp=(1:fix(J/2)+1)/(J*dt);
X=plomb(x-mxe,t,fp)/2;
E=zeros(1,J);
for k=2:fix(J/2)+1
E(k)=T^2*(N^2*X(k-1)/T)/(N^2);
end
for k=fix(J/2)+2:J
E(k)=E(J-k+2);
end
RE=real(ifft(E))/dt;
K=200;
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros(1,K);
for k=1:fix(K/2)+1
R(k)=RE(k)/(T-abs(tau(k)))+mxe^2;
end
for k=fix(K/2)+2:K
R(k)=R(K-k+2);
end
plot(tau,R,',',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs^2)
pause
%
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'.',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über direkte Spektralschätzung (Lomb-Scargle-Methode, mit Abzug der Selbstprodukte, ohne Normierung, ohne Korrektur des dynamischen Fehlers)
from numpy.fft import *
import scipy.signal as scisig
N=len(x)
dt=1.0
mxe=mean(x)
J=int(ceil(2*T/float(dt)))
fp=arange(1,J//2+1)/float(J*dt)
X=scisig.lombscargle(t,x-mxe,2*pi*fp)
E=zeros(J)
for k in range(0,J//2+1):
  E[k]=T**2*(N*X[k-1]-sum((x-mxe)**2))/float(N**2-N)

for k in range(-((J-1)//2),0):
  E[k]=E[-k]

RE=real(ifft(E))/float(dt)
K=200
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)
for k in range(-(K//2),(K+1)//2):
  R[k]=RE[k]/float(T-abs(tau[k]))+mxe**2

plot(tau,R,'.',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs**2)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'.',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(x);
dt=1.0;
mxe=mean(x);
J=ceil(2*T/dt);
fp=(1:fix(J/2)+1)/(J*dt);
X=plomb(x-mxe,t,fp)/2;
E=zeros(1,J);
for k=2:fix(J/2)+1
E(k)=T^2*(N^2*X(k-1)/T-sum((x-mxe).^2))/(N^2-N);
end
for k=fix(J/2)+2:J
E(k)=E(J-k+2);
end
RE=real(ifft(E))/dt;
K=200;
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros(1,K);
for k=1:fix(K/2)+1
R(k)=RE(k)/(T-abs(tau(k)))+mxe^2;
end
for k=fix(K/2)+2:K
R(k)=R(K-k+2);
end
plot(tau,R,'.',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs^2)
pause
%
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'.',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Zeitquantisierung (mit Abzug der Selbstprodukte, mit Normierung)
N=len(t)
dt=1.0
J=int(ceil(2*T/float(dt)))
x1=zeros(J)
x0=zeros(J)
for i in range(0,N):
  j=int(floor(t[i]/float(dt)))
  if j<J:
    x1[j]+=x[i];
    x0[j]+=1;
  

X=fft(x1)
Xp=fft(x0)
E=T**2*(abs(X)**2-sum(x**2))/float(N**2-N)
Ep=T**2*(abs(Xp)**2-N)/float(N**2-N)
RE=real(ifft(E))/float(dt)
REp=real(ifft(Ep))/float(dt)
K=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)
for k in range(-(K//2),(K+1)//2):
  if REp[k]!=0:
    R[k]=RE[k]/REp[k]
  

plot(tau,R,'.',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs**2)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'.',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
J=ceil(2*T/dt);
x1=zeros(1,J);
x0=zeros(1,J);
for i=1:N
j=fix(t(i)/dt)+1;
if j<=J
x1(j)=x1(j)+x(i);
x0(j)=x0(j)+1;
end
end
X=fft(x1);
Xp=fft(x0);
E=T^2*(abs(X).^2-sum(x.^2))/(N^2-N);
Ep=T^2*(abs(Xp).^2-N)/(N^2-N);
RE=real(ifft(E))/dt;
REp=real(ifft(Ep))/dt;
K=200;
K=min(ceil(2*T/dt)-1,K);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros(1,K);
for k=1:fix(K/2)+1
if REp(k)~=0
R(k)=RE(k)/REp(k);
end
end
for k=fix(K/2)+2:K
R(k)=R(K-k+2);
end
plot(tau,R,'.',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs^2)
pause
%
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'.',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Interpolation (Sample-and-Hold, ohne Filterkorrektur, ohne Normierung)
N=len(t)
dt=1.0
J=int(ceil(2*T/float(dt)))
xp=zeros(J)
for i in range(0,N):
  for j in range(int(floor(t[i]/float(dt)))+1,int(floor((t[(i+1)%N]+T*((i+1)//N))/float(dt)))+1):
    xp[j]=x[i]
  

Xp=fft(xp)
Ep=dt**2*abs(Xp)**2
REp=real(ifft(Ep))/float(dt)
K=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)
for k in range(-(K//2),(K+1)//2):
  R[k]=REp[k]/float(T-abs(tau[k]))

plot(tau,R,'.',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs**2)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'.',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
J=ceil(2*T/dt);
xp=zeros(1,J);
for i=1:N
for j=fix(t(i)/dt)+2:fix((t(mod(i,N)+1)+T*fix(i/N))/dt)+1
xp(j)=x(i);
end
end
Xp=fft(xp);
Ep=dt^2*abs(Xp).^2;
REp=real(ifft(Ep))/dt;
K=200;
K=min(ceil(2*T/dt)-1,K);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros(1,K);
for k=1:fix(K/2)+1
R(k)=REp(k)/(T-abs(tau(k)));
end
for k=fix(K/2)+2:K
R(k)=R(K-k+2);
end
plot(tau,R,'.',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs^2)
pause
%
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'.',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Interpolation (Sample-and-Hold, mit Filterkorrektur, ohne Normierung)
N=len(t)
dt=1.0
J=int(ceil(2*T/float(dt)))
xp=zeros(J)
for i in range(0,N):
  for j in range(int(floor(t[i]/float(dt)))+1,int(floor((t[(i+1)%N]+T*((i+1)//N))/float(dt)))+1):
    xp[j]=x[i]
  

Xp=fft(xp)
Ep=dt**2*abs(Xp)**2
REp=real(ifft(Ep))/float(dt)
K=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)
c=exp(-dd/float(dt))/(1-exp(-dd/float(dt)))**2
for k in range(-(K//2),(K+1)//2):
  if k==0:
    R[k]=REp[k]/float(T)
  else:
    R[k]=((2*c+1)*REp[k]-c*REp[k-1]-c*REp[k+1])/float(T-abs(tau[k]))
  

plot(tau,R,'.',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs**2)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'.',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
J=ceil(2*T/dt);
xp=zeros(1,J);
for i=1:N
for j=fix(t(i)/dt)+2:fix((t(mod(i,N)+1)+T*fix(i/N))/dt)+1
xp(j)=x(i);
end
end
Xp=fft(xp);
Ep=dt^2*abs(Xp).^2;
REp=real(ifft(Ep))/dt;
K=200;
K=min(ceil(2*T/dt)-1,K);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros(1,K);
c=exp(-dd/dt)/(1-exp(-dd/dt))^2;
R(1)=REp(1)/T;
for k=2:fix(K/2)+1
R(k)=((2*c+1)*REp(k)-c*REp(k-1)-c*REp(k+1))/(T-abs(tau(k)));
end
for k=fix(K/2)+2:K
R(k)=R(K-k+2);
end
plot(tau,R,'.',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs^2)
pause
%
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'.',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Generieren von Wertepaaren (ti,xi) mit den Messzeitpunkten ti und den Messwerten xi (AR1-Prozess als Beispiel, mit Mittelwert verschieden von null, neuer Datensatz mit geringerer Datenrate)
T=10000.0
Ti=10.0
dd=0.1
mxs=3.0
vxs=1.0
t=[]
x=[]
te=0.0
xe=normal(mxs,sqrt(vxs))
while te<T:
  tp=exponential(1.0/dd)
  te+=tp
  if te<T:
    t.append(te)
    phi=exp(-tp/Ti)
    theta=sqrt((1-phi**2)*vxs)
    xe=(xe-mxs)*phi+normal(mxs,theta)
    x.append(xe)
  

t=array(t)
x=array(x)
plot(t,x,'.')
show()
T=10000;
Ti=10;
dd=0.1;
mxs=3;
vxs=1;
t=[];
x=[];
te=0;
xe=normrnd(mxs,sqrt(vxs));
while te<T
tp=exprnd(1/dd);
te=te+tp;
if te<T
t(end+1)=te;
phi=exp(-tp/Ti);
theta=sqrt((1-phi^2)*vxs);
xe=(xe-mxs)*phi+normrnd(mxs,theta);
x(end+1)=xe;
end
end
plot(t,x,'.')
Autokorrelationsfunktion und Leistungsdichtespektrum über Interpolation (Sample-and-Hold, ohne Filterkorrektur, ohne Normierung)
N=len(t)
dt=1.0
J=int(ceil(2*T/float(dt)))
xp=zeros(J)
for i in range(0,N):
  for j in range(int(floor(t[i]/float(dt)))+1,int(floor((t[(i+1)%N]+T*((i+1)//N))/float(dt)))+1):
    xp[j]=x[i]
  

Xp=fft(xp)
Ep=dt**2*abs(Xp)**2
REp=real(ifft(Ep))/float(dt)
K=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)
for k in range(-(K//2),(K+1)//2):
  R[k]=REp[k]/float(T-abs(tau[k]))

plot(tau,R,'.',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs**2)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'.',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
J=ceil(2*T/dt);
xp=zeros(1,J);
for i=1:N
for j=fix(t(i)/dt)+2:fix((t(mod(i,N)+1)+T*fix(i/N))/dt)+1
xp(j)=x(i);
end
end
Xp=fft(xp);
Ep=dt^2*abs(Xp).^2;
REp=real(ifft(Ep))/dt;
K=200;
K=min(ceil(2*T/dt)-1,K);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros(1,K);
for k=1:fix(K/2)+1
R(k)=REp(k)/(T-abs(tau(k)));
end
for k=fix(K/2)+2:K
R(k)=R(K-k+2);
end
plot(tau,R,'.',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs^2)
pause
%
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'.',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))
Autokorrelationsfunktion und Leistungsdichtespektrum über Interpolation (Sample-and-Hold, mit Filterkorrektur, ohne Normierung)
N=len(t)
dt=1.0
J=int(ceil(2*T/float(dt)))
xp=zeros(J)
for i in range(0,N):
  for j in range(int(floor(t[i]/float(dt)))+1,int(floor((t[(i+1)%N]+T*((i+1)//N))/float(dt)))+1):
    xp[j]=x[i]
  

Xp=fft(xp)
Ep=dt**2*abs(Xp)**2
REp=real(ifft(Ep))/float(dt)
K=200
K=minimum(int(ceil(2*T/float(dt)))-1,K)
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
R=zeros(K)
c=exp(-dd/float(dt))/(1-exp(-dd/float(dt)))**2
for k in range(-(K//2),(K+1)//2):
  if k==0:
    R[k]=REp[k]/float(T)
  else:
    R[k]=((2*c+1)*REp[k]-c*REp[k-1]-c*REp[k+1])/float(T-abs(tau[k]))
  

plot(tau,R,'.',arange(-5*K,5*K)*dt/10.0,vxs*exp(-abs(arange(-5*K,5*K)*dt/10.0/float(Ti)))+mxs**2)
show()
f=roll(arange(-(K//2),(K+1)//2)/float(K*dt),(K+1)//2)
S=dt*real(fft(R))
p1=1-dt/float(Ti)
loglog(f,S,'.',arange(1,5*K)/float(10*K*dt),vxs*(1-p1**2)*dt/(1+p1**2-2*p1*cos(2*pi*arange(1,5*K)/float(10*K*dt)*dt)))
show()
N=length(t);
dt=1.0;
J=ceil(2*T/dt);
xp=zeros(1,J);
for i=1:N
for j=fix(t(i)/dt)+2:fix((t(mod(i,N)+1)+T*fix(i/N))/dt)+1
xp(j)=x(i);
end
end
Xp=fft(xp);
Ep=dt^2*abs(Xp).^2;
REp=real(ifft(Ep))/dt;
K=200;
K=min(ceil(2*T/dt)-1,K);
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,[0;fix((K+1)/2)]);
R=zeros(1,K);
c=exp(-dd/dt)/(1-exp(-dd/dt))^2;
R(1)=REp(1)/T;
for k=2:fix(K/2)+1
R(k)=((2*c+1)*REp(k)-c*REp(k-1)-c*REp(k+1))/(T-abs(tau(k)));
end
for k=fix(K/2)+2:K
R(k)=R(K-k+2);
end
plot(tau,R,'.',(-5*K:5*K)*dt/10,vxs*exp(-abs((-5*K:5*K)*dt/10/Ti))+mxs^2)
pause
%
f=circshift((-fix(K/2):fix((K-1)/2))/(K*dt),[0;fix((K+1)/2)]);
S=dt*real(fft(R));
p1=1-dt/Ti;
loglog(f,S,'.',(1:5*K)/(10*K*dt),vxs*(1-p1^2)*dt./(1+p1^2-2*p1*cos(2*pi*(1:5*K)/(10*K*dt)*dt)))