Signal- und Messdatenverarbeitung

Codes für komplexe, 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 komplexen Werten xi mit i=0…N−1 (harmonisches Signal als Beispiel, mit einem Erwartungswert verschieden von null, ohne Rauschen)
N=100
dt=1.0
t=arange(0,N)*dt
x=3+1j+1.5*exp(0.1j*pi*t)
plot(t,real(x),'o',t,imag(x),'o')
show()
N=100;
dt=1.0;
t=(0:N-1)*dt;
x=3+1j+1.5*exp(0.1j*pi*t);
plot(t,real(x),'o',t,imag(x),'o')
Überlagerung von Rauschen
x+=normal(0,0.1,N)+1j*normal(0,0.1,N)
plot(t,real(x),'o',t,imag(x),'o')
show()
x=x+0.1*randn(1,N)+1j*0.1*randn(1,N);
plot(t,real(x),'o',t,imag(x),'o')
Mittelwert (ohne Rauschen exakt, mit Rauschen erwartungstreu) x‾=1N⋅∑i=0N−1xi
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=1N⋅∑i=0N−1(xi−x‾)*(xi−x‾)=1N⋅∑i=0N−1|xi−x‾|2
=(1N⋅∑i=0N−1xi*xi)−x‾*x‾=(1N⋅∑i=0N−1|xi|2)−|x‾|2
var(x)
var(x,1)
(Auto-)Korrelationsfunktion (ohne Rauschen exakt, mit Rauschen systematischer Fehler bei τk=0) Rk=R(τk)=1N⋅∑i=0N−1xi*x(i+k)modN
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
R=zeros(N)+0j
tau=arange(0,N)*dt
for k in range(0,N):
  sum1=0+0j
  for i in range(0,N):
    sum1+=conj(x[i])*x[(i+k)%N]
  
  R[k]=sum1/N

plot(tau,real(R),'o',tau,imag(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*ifft(P)
plot(tau,real(R),'o',tau,imag(R),'o')
show()
 direkte Berechnung:
R=zeros(1,N)+0j;
tau=(0:N-1)*dt;
for k=0:N-1
sum1=0+0j;
for i=1:N
sum1=sum1+conj(x(i))*x(mod(i+k-1,N)+1);
end
R(k+1)=sum1/N;
end
plot(tau,real(R),'o',tau,imag(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*ifft(P);
plot(tau,real(R),'o',tau,imag(R),'o')
Berechnung für den Fall, dass die anzuwendende Periode des Signals Ñ von der Signallänge N abweicht (Ñ≤N) Rk=R(τk)=∑i=0N−1∑j=0N−1j−i≡kmodÑxi*xj∑i=0N−1∑j=0N−1j−i≡kmodÑ1
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
NT=60 #Periodenlänge NT≤N
R1=zeros(NT)+0j
R0=zeros(NT)
for i in range(0,N):
  for j in range(0,N):
    R1[(j-i)%NT]+=conj(x[i])*x[j]
    R0[(j-i)%NT]+=1
  

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

plot(tau,real(R),'o',tau,imag(R),'o')
show()
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
NT=60 #Periodenlänge NT≤N
xp=zeros(NT)+0j
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=ifft(P1)/real(ifft(P0))
plot(tau,real(R),'o',tau,imag(R),'o')
show()
 direkte Berechnung:
NT=60; %Periodenlänge NT≤N
R1=zeros(1,NT)+0j;
R0=zeros(1,NT);
for i=1:N
for j=1:N
R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+conj(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)+0j;
for k=1:NT
R(k)=R1(k)/R0(k);
end
plot(tau,real(R),'o',tau,imag(R),'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
NT=60; %Periodenlänge NT≤N
xp=zeros(1,NT)+0j;
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=ifft(P1)./real(ifft(P0));
plot(tau,real(R),'o',tau,imag(R),'o')
(Auto-)Kovarianzfunktion (ohne Rauschen exakt, mit Rauschen systematischer Fehler bei τk=0 und bei τk≠0 asymptotisch erwatungstreu) Ck=C(τk)=1N⋅∑i=0N−1(xi−x‾)*[x(i+k)modN−x‾]
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
C=zeros(N)+0j
tau=arange(0,N)*dt
mxe=mean(x)
for k in range(0,N):
  sum1=0+0j
  for i in range(0,N):
    sum1+=conj(x[i]-mxe)*(x[(i+k)%N]-mxe)
  
  C[k]=sum1/N

plot(tau,real(C),'o',tau,imag(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*ifft(P)
plot(tau,real(C),'o',tau,imag(C),'o')
show()
 direkte Berechnung:
C=zeros(1,N)+0j;
tau=(0:N-1)*dt;
mxe=mean(x);
for k=0:N-1
sum1=0+0j;
for i=1:N
sum1=sum1+conj(x(i)-mxe)*(x(mod(i+k-1,N)+1)-mxe);
end
C(k+1)=sum1/N;
end
plot(tau,real(C),'o',tau,imag(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*ifft(P);
plot(tau,real(C),'o',tau,imag(C),'o')
Berechnung für den Fall, dass die anzuwendende Periode des Signals Ñ von der Signallänge N abweicht (Ñ≤N) Ck=C(τk)=∑i=0N−1∑j=0N−1j−i≡kmodÑ(xi−x‾)*(xj−x‾)∑i=0N−1∑j=0N−1j−i≡kmodÑ1
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
mxe=mean(x)
NT=60 #Periodenlänge NT≤N
C1=zeros(NT)+0j
C0=zeros(NT)
for i in range(0,N):
  for j in range(0,N):
    C1[(j-i)%NT]+=conj(x[i]-mxe)*(x[j]-mxe)
    C0[(j-i)%NT]+=1
  

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

plot(tau,real(C),'o',tau,imag(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)+0j
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=ifft(P1)/real(ifft(P0))
plot(tau,real(C),'o',tau,imag(C),'o')
show()
 direkte Berechnung:
mxe=mean(x);
NT=60; %Periodenlänge NT≤N
C1=zeros(1,NT)+0j;
C0=zeros(1,NT);
for i=1:N
for j=1:N
C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+conj(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)+0j;
for k=1:NT
C(k)=C1(k)/C0(k);
end
plot(tau,real(C),'o',tau,imag(C),'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
mxe=mean(x);
NT=60; %Periodenlänge NT≤N
xp=zeros(1,NT)+0j;
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=ifft(P1)./real(ifft(P0));
plot(tau,real(C),'o',tau,imag(C),'o')
(Auto-)Korrelationskoeffizientenfunktion (ohne Rauschen exakt, mit Rauschen systematischer Fehler bei allen τk≠0) ρk=ρ(τk)=Cks2=1Ns2⋅∑i=0N−1(xi−x‾)*[x(i+k)modN−x‾]
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
rho=zeros(N)+0j
tau=arange(0,N)*dt
mxe=mean(x)
vxe=var(x)
for k in range(0,N):
  sum1=0+0j
  for i in range(0,N):
    sum1+=conj(x[i]-mxe)*(x[(i+k)%N]-mxe)
  
  rho[k]=sum1/N/vxe

plot(tau,real(rho),'o',tau,imag(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*ifft(P)
rho=C/var(x)
plot(tau,real(rho),'o',tau,imag(rho),'o')
show()
 direkte Berechnung:
rho=zeros(1,N)+0j;
tau=(0:N-1)*dt;
mxe=mean(x);
vxe=var(x,1);
for k=0:N-1
sum1=0+0j;
for i=1:N
sum1=sum1+conj(x(i)-mxe)*(x(mod(i+k-1,N)+1)-mxe);
end
rho(k+1)=sum1/N/vxe;
end
plot(tau,real(rho),'o',tau,imag(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*ifft(P);
rho=C/var(x,1);
plot(tau,real(rho),'o',tau,imag(rho),'o')
Berechnung für den Fall, dass die anzuwendende Periode des Signals Ñ von der Signallänge N abweicht (Ñ≤N) ρk=ρ(τk)=Cks2=∑i=0N−1∑j=0N−1j−i≡kmodÑ(xi−x‾)*(xj−x‾)s2⋅∑i=0N−1∑j=0N−1j−i≡kmodÑ1
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
mxe=mean(x)
vxe=var(x)
NT=60 #Periodenlänge NT≤N
C1=zeros(NT)+0j
C0=zeros(NT)
for i in range(0,N):
  for j in range(0,N):
    C1[(j-i)%NT]+=conj(x[i]-mxe)*(x[j]-mxe)
    C0[(j-i)%NT]+=1
  

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

plot(tau,real(rho),'o',tau,imag(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)+0j
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=ifft(P1)/real(ifft(P0))/vxe
plot(tau,real(rho),'o',tau,imag(rho),'o')
show()
 direkte Berechnung:
mxe=mean(x);
vxe=var(x,1);
NT=60; %Periodenlänge NT≤N
C1=zeros(1,NT)+0j;
C0=zeros(1,NT);
for i=1:N
for j=1:N
C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+conj(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)+0j;
for k=1:NT
rho(k)=C1(k)/C0(k)/vxe;
end
plot(tau,real(rho),'o',tau,imag(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)+0j;
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=ifft(P1)./real(ifft(P0))/vxe;
plot(tau,real(rho),'o',tau,imag(rho),'o')
(Auto-)Leistungsspektrum (ohne Rauschen exakt, mit Rauschen systematischer Fehler bei allen Frequenzen) Xj=FFT{xi}=∑i=0N−1xi𝐞−2π𝐢ij/N
(imaginäre Einheit 𝐢)
Pj=P(fj)=Xj*XjN2=|Xj|2N2
fj=jNΔt
j=−⌊N/2⌋…⌊(N−1)/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)+0j
for k in range(0,N):
  sum1=0+0j
  for i in range(0,N):
    sum1+=conj(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)+0j;
for k=0:N-1
sum1=0+0j;
for i=1:N
sum1=sum1+conj(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)+0j
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=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)+0j
R0=zeros(NT)
for i in range(0,N):
  for j in range(0,N):
    R1[(j-i)%NT]+=conj(x[i])*x[j]
    R0[(j-i)%NT]+=1
  

R=zeros(NT)+0j
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)+0j;
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=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)+0j;
R0=zeros(1,NT);
for i=1:N
for j=1:N
R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+conj(x(i))*x(j);
R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+1;
end
end
R=zeros(1,NT)+0j;
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')