Signal- und Messdatenverarbeitung


Codes für reelle, periodische Einzelsignale, ohne Gewichtung

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
Generieren von N Werten xi mit i=0N1 (Superposition zweier harmonischer Signale als Beispiel, mit einem Erwartungswert verschieden von null, ohne Rauschen)
N=100
dt=1.0
t=arange(0,N)*dt
x=3+1.5*sin(0.2*pi*t)+1.5*sin(0.1*pi*t)
plot(t,x,'o')
show()
N=100;
dt=1.0;
t=(0:N-1)*dt;
x=3+1.5*sin(0.2*pi*t)+1.5*sin(0.1*pi*t);
plot(t,x,'o')
Überlagerung von Rauschen
x+=normal(0,0.1,N)
plot(t,x,'o')
show()
x=x+0.1*randn(1,N);
plot(t,x,'o')
Mittelwert (ohne Rauschen exakt, mit Rauschen erwartungstreu) x=1Ni=0N1xi
mean(x)
mean(x)
Varianz (für streng periodische Signale keine Bessel-Korrektur, da Mittelwert keine Schätzgröße, ohne Rauschen exakt, mit Rauschen systematischer Fehler) s2=1Ni=0N1(xix)2=(1Ni=0N1xi2)x2
var(x)
var(x,1)
(Auto-)Korrelationsfunktion (ohne Rauschen exakt, mit Rauschen systematischer Fehler bei τk=0) Rk=R(τk)=1Ni=0N1xix(i+k)modN
τk=kΔt
k=
direkte Berechnung:
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/N

plot(tau,R,'o')
show()
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
tau=arange(0,N)*dt
X=fft(x)
P=abs(X)**2/N**2
R=N*real(ifft(P))
plot(tau,R,'o')
show()
direkte Berechnung:
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,'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
tau=(0:N-1)*dt;
X=fft(x);
P=abs(X).^2/N^2;
R=N*real(ifft(P));
plot(tau,R,'o')
Berechnung für den Fall, dass die anzuwendende Periode des Signals Ñ von der Signallänge N abweicht (ÑN) Rk=R(τk)=i=0N1j=0N1jikmodÑxixji=0N1j=0N1jikmodÑ1
τk=kΔt
k=
direkte Berechnung:
NT=60 #Periodenlänge NT≤N
R1=zeros(NT)
R0=zeros(NT)
for i in range(0,N):
  for j in range(0,N):
    R1[(j-i)%NT]+=x[i]*x[j]
    R0[(j-i)%NT]+=1
  

tau=arange(0,NT)*dt
R=zeros(NT)
for k in range(0,NT):
  R[k]=R1[k]/R0[k]

plot(tau,R,'o')
show()
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
NT=60 #Periodenlänge NT≤N
xp=zeros(NT)
np=zeros(NT)
for i in range(0,N):
  xp[i%NT]+=x[i]
  np[i%NT]+=1

tau=arange(0,NT)*dt
X1=fft(xp)
X0=fft(np)
P1=abs(X1)**2/NT**2
P0=abs(X0)**2/NT**2
R=real(ifft(P1))/real(ifft(P0))
plot(tau,R,'o')
show()
direkte Berechnung:
NT=60; %Periodenlänge NT≤N
R1=zeros(1,NT);
R0=zeros(1,NT);
for i=1:N
for j=1:N
R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+x(i)*x(j);
R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+1;
end
end
tau=(0:NT-1)*dt;
R=zeros(1,NT);
for k=1:NT
R(k)=R1(k)/R0(k);
end
plot(tau,R,'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
NT=60; %Periodenlänge NT≤N
xp=zeros(1,NT);
np=zeros(1,NT);
for i=1:N
xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+x(i);
np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1;
end
tau=(0:NT-1)*dt;
X1=fft(xp);
X0=fft(np);
P1=abs(X1).^2/NT^2;
P0=abs(X0).^2/NT^2;
R=real(ifft(P1))./real(ifft(P0));
plot(tau,R,'o')
(Auto-)Kovarianzfunktion (ohne Rauschen exakt, mit Rauschen systematischer Fehler bei τk=0 und bei τk0 asymptotisch erwatungstreu) Ck=C(τk)=1Ni=0N1(xix)[x(i+k)modNx]
τk=kΔt
k=
direkte Berechnung:
C=zeros(N)
tau=arange(0,N)*dt
mxe=mean(x)
for k in range(0,N):
  sum1=0
  for i in range(0,N):
    sum1+=(x[i]-mxe)*(x[(i+k)%N]-mxe)
  
  C[k]=sum1/N

plot(tau,C,'o')
show()
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
tau=arange(0,N)*dt
X=fft(x)
P=abs(X)**2/N**2
P[0]=0
C=N*real(ifft(P))
plot(tau,C,'o')
show()
direkte Berechnung:
C=zeros(1,N);
tau=(0:N-1)*dt;
mxe=mean(x);
for k=0:N-1
sum1=0;
for i=1:N
sum1=sum1+(x(i)-mxe)*(x(mod(i+k-1,N)+1)-mxe);
end
C(k+1)=sum1/N;
end
plot(tau,C,'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
tau=(0:N-1)*dt;
X=fft(x);
P=abs(X).^2/N^2;
P(1)=0;
C=N*real(ifft(P));
plot(tau,C,'o')
Berechnung für den Fall, dass die anzuwendende Periode des Signals Ñ von der Signallänge N abweicht (ÑN) Ck=C(τk)=i=0N1j=0N1jikmodÑ(xix)(xjx)i=0N1j=0N1jikmodÑ1
τk=kΔt
k=
direkte Berechnung:
mxe=mean(x)
NT=60 #Periodenlänge NT≤N
C1=zeros(NT)
C0=zeros(NT)
for i in range(0,N):
  for j in range(0,N):
    C1[(j-i)%NT]+=(x[i]-mxe)*(x[j]-mxe)
    C0[(j-i)%NT]+=1
  

tau=arange(0,NT)*dt
C=zeros(NT)
for k in range(0,NT):
  C[k]=C1[k]/C0[k]

plot(tau,C,'o')
show()
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
mxe=mean(x)
NT=60 #Periodenlänge NT≤N
xp=zeros(NT)
np=zeros(NT)
for i in range(0,N):
  xp[i%NT]+=x[i]-mxe
  np[i%NT]+=1

tau=arange(0,NT)*dt
X1=fft(xp)
X0=fft(np)
P1=abs(X1)**2/NT**2
P0=abs(X0)**2/NT**2
C=real(ifft(P1))/real(ifft(P0))
plot(tau,C,'o')
show()
direkte Berechnung:
mxe=mean(x);
NT=60; %Periodenlänge NT≤N
C1=zeros(1,NT);
C0=zeros(1,NT);
for i=1:N
for j=1:N
C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+(x(i)-mxe)*(x(j)-mxe);
C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+1;
end
end
tau=(0:NT-1)*dt;
C=zeros(1,NT);
for k=1:NT
C(k)=C1(k)/C0(k);
end
plot(tau,C,'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
mxe=mean(x);
NT=60; %Periodenlänge NT≤N
xp=zeros(1,NT);
np=zeros(1,NT);
for i=1:N
xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+x(i)-mxe;
np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1;
end
tau=(0:NT-1)*dt;
X1=fft(xp);
X0=fft(np);
P1=abs(X1).^2/NT^2;
P0=abs(X0).^2/NT^2;
C=real(ifft(P1))./real(ifft(P0));
plot(tau,C,'o')
(Auto-)Korrelationskoeffizientenfunktion (ohne Rauschen exakt, mit Rauschen systematischer Fehler bei allen τk0) ρk=ρ(τk)=Cks2=1Ns2i=0N1(xix)[x(i+k)modNx]
τk=kΔt
k=
direkte Berechnung:
rho=zeros(N)
tau=arange(0,N)*dt
mxe=mean(x)
vxe=var(x)
for k in range(0,N):
  sum1=0
  for i in range(0,N):
    sum1+=(x[i]-mxe)*(x[(i+k)%N]-mxe)
  
  rho[k]=sum1/N/vxe

plot(tau,rho,'o')
show()
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
tau=arange(0,N)*dt
X=fft(x)
P=abs(X)**2/N**2
P[0]=0
C=N*real(ifft(P))
rho=C/var(x)
plot(tau,rho,'o')
show()
direkte Berechnung:
rho=zeros(1,N);
tau=(0:N-1)*dt;
mxe=mean(x);
vxe=var(x,1);
for k=0:N-1
sum1=0;
for i=1:N
sum1=sum1+(x(i)-mxe)*(x(mod(i+k-1,N)+1)-mxe);
end
rho(k+1)=sum1/N/vxe;
end
plot(tau,rho,'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
tau=(0:N-1)*dt;
X=fft(x);
P=abs(X).^2/N^2;
P(1)=0;
C=N*real(ifft(P));
rho=C/var(x,1);
plot(tau,rho,'o')
Berechnung für den Fall, dass die anzuwendende Periode des Signals Ñ von der Signallänge N abweicht (ÑN) ρk=ρ(τk)=Cks2=i=0N1j=0N1jikmodÑ(xix)(xjx)s2i=0N1j=0N1jikmodÑ1
τk=kΔt
k=
direkte Berechnung:
mxe=mean(x)
vxe=var(x)
NT=60 #Periodenlänge NT≤N
C1=zeros(NT)
C0=zeros(NT)
for i in range(0,N):
  for j in range(0,N):
    C1[(j-i)%NT]+=(x[i]-mxe)*(x[j]-mxe)
    C0[(j-i)%NT]+=1
  

tau=arange(0,NT)*dt
rho=zeros(NT)
for k in range(0,NT):
  rho[k]=C1[k]/C0[k]/vxe

plot(tau,rho,'o')
show()
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
mxe=mean(x)
vxe=var(x)
NT=60 #Periodenlänge NT≤N
xp=zeros(NT)
np=zeros(NT)
for i in range(0,N):
  xp[i%NT]+=x[i]-mxe
  np[i%NT]+=1

tau=arange(0,NT)*dt
X1=fft(xp)
X0=fft(np)
P1=abs(X1)**2/NT**2
P0=abs(X0)**2/NT**2
rho=real(ifft(P1))/real(ifft(P0))/vxe
plot(tau,rho,'o')
show()
direkte Berechnung:
mxe=mean(x);
vxe=var(x,1);
NT=60; %Periodenlänge NT≤N
C1=zeros(1,NT);
C0=zeros(1,NT);
for i=1:N
for j=1:N
C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+(x(i)-mxe)*(x(j)-mxe);
C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+1;
end
end
tau=(0:NT-1)*dt;
rho=zeros(1,NT);
for k=1:NT
rho(k)=C1(k)/C0(k)/vxe;
end
plot(tau,rho,'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
mxe=mean(x);
vxe=var(x,1);
NT=60; %Periodenlänge NT≤N
xp=zeros(1,NT);
np=zeros(1,NT);
for i=1:N
xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+x(i)-mxe;
np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1;
end
tau=(0:NT-1)*dt;
X1=fft(xp);
X0=fft(np);
P1=abs(X1).^2/NT^2;
P0=abs(X0).^2/NT^2;
rho=real(ifft(P1))./real(ifft(P0))/vxe;
plot(tau,rho,'o')
(Auto-)Leistungsspektrum (ohne Rauschen exakt, mit Rauschen systematischer Fehler bei allen Frequenzen) Xj=FFT{xi}=i=0N1xi𝐞2π𝐢ij/N
(imaginäre Einheit 𝐢)
Pj=P(fj)=|Xj|2N2
fj=jNΔt
j=N/2(N1)/2
direkte Berechnung:
from numpy.fft import *
f=roll(arange(-(N//2),(N+1)//2)/float(N*dt),(N+1)//2)
X=fft(x)
P=abs(X)**2/N**2
plot(f,P,'o')
show()
Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
R=zeros(N)
for k in range(0,N):
  sum1=0
  for i in range(0,N):
    sum1+=x[i]*x[(i+k)%N]
  
  R[k]=sum1/N

f=roll(arange(-(N//2),(N+1)//2)/float(N*dt),(N+1)//2)
P=real(fft(R))/float(N)
plot(f,P,'o')
show()
direkte Berechnung:
f=circshift((-fix(N/2):fix((N-1)/2))/(N*dt),[0;fix((N+1)/2)]);
X=fft(x);
P=abs(X).^2/N^2;
plot(f,P,'o')
Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
R=zeros(1,N);
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
f=circshift((-fix(N/2):fix((N-1)/2))/(N*dt),[0;fix((N+1)/2)]);
P=real(fft(R))/N;
plot(f,P,'o')
Berechnung für den Fall, dass die anzuwendende Periode des Signals Ñ von der Signallänge N abweicht (ÑN) (wegen der nötigen Entfaltung mit dem Leistungsspektrum der Überlappungen nur über die Korrelationsfunktion)
Pj=P(fj)=1ÑFFT{Rk}=1Ñk=0Ñ1Rk𝐞2π𝐢jk/Ñ
(imaginäre Einheit 𝐢)
fj=jÑΔt
j=Ñ/2(Ñ1)/2
Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren:
from numpy.fft import *
NT=60 #Periodenlänge NT≤N
xp=zeros(NT)
np=zeros(NT)
for i in range(0,N):
  xp[i%NT]+=x[i]
  np[i%NT]+=1

X1=fft(xp)
X0=fft(np)
P1=abs(X1)**2/NT**2
P0=abs(X0)**2/NT**2
R=real(ifft(P1))/real(ifft(P0))
f=roll(arange(-(NT//2),(NT+1)//2)/float(NT*dt),(NT+1)//2)
P=real(fft(R))/float(NT)
plot(f,P,'o')
show()
Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
NT=60 #Periodenlänge NT≤N
R1=zeros(NT)
R0=zeros(NT)
for i in range(0,N):
  for j in range(0,N):
    R1[(j-i)%NT]+=x[i]*x[j]
    R0[(j-i)%NT]+=1
  

R=zeros(NT)
for k in range(0,NT):
  R[k]=R1[k]/R0[k]

f=roll(arange(-(NT//2),(NT+1)//2)/float(NT*dt),(NT+1)//2)
P=real(fft(R))/float(NT)
plot(f,P,'o')
show()
Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren:
NT=60; %Periodenlänge NT≤N
xp=zeros(1,NT);
np=zeros(1,NT);
for i=1:N
xp(mod(i-1,NT)+1)=xp(mod(i-1,NT)+1)+x(i);
np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1;
end
X1=fft(xp);
X0=fft(np);
P1=abs(X1).^2/NT^2;
P0=abs(X0).^2/NT^2;
R=real(ifft(P1))./real(ifft(P0));
f=circshift((-fix(NT/2):fix((NT-1)/2))/(NT*dt),[0;fix((NT+1)/2)]);
P=real(fft(R))/NT;
plot(f,P,'o')
Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
NT=60; %Periodenlänge NT≤N
R1=zeros(1,NT);
R0=zeros(1,NT);
for i=1:N
for j=1:N
R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+x(i)*x(j);
R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+1;
end
end
R=zeros(1,NT);
for k=1:NT
R(k)=R1(k)/R0(k);
end
f=circshift((-fix(NT/2):fix((NT-1)/2))/(NT*dt),[0;fix((NT+1)/2)]);
P=real(fft(R))/NT;
plot(f,P,'o')