Signal- und Messdatenverarbeitung

Codes für komplexe, periodische Signalpaare, 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 Wertepaaren (xi,yi) 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)
y=2+1.5j+1*exp(0.1j*pi*t-0.3)
plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o')
show()
N=100;
dt=1.0;
t=(0:N-1)*dt;
x=3+1j+1.5*exp(0.1j*pi*t);
y=2+1.5j+1*exp(0.1j*pi*t-0.3);
plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o')
Überlagerung von Rauschen
x+=normal(0,0.1,N)+1j*normal(0,0.1,N)
y+=normal(0,0.1,N)+1j*normal(0,0.1,N)
plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o')
show()
x=x+0.1*randn(1,N)+1j*0.1*randn(1,N);
y=y+0.1*randn(1,N)+1j*0.1*randn(1,N);
plot(t,real(x),'o',t,imag(x),'o',t,real(y),'o',t,imag(y),'o')
Mittelwerte (ohne Rauschen exakt, mit Rauschen erwartungstreu) x‾=1N⋅∑i=0N−1xi
y‾=1N⋅∑i=0N−1yi
mxe=mean(x)
mye=mean(y)
mxe=mean(x);
mye=mean(y);
Varianzen (für streng periodische Signale keine Bessel-Korrektur, da Mittelwert keine Schätzgröße, ohne Rauschen exakt, mit Rauschen systematischer Fehler) sx2=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
sy2=1N⋅∑i=0N−1(yi−y‾)*(yi−y‾)=1N⋅∑i=0N−1|yi−y‾|2
=(1N⋅∑i=0N−1yi*yi)−y‾*y‾=(1N⋅∑i=0N−1|yi|2)−|y‾|2
vxe=var(x)
vye=var(y)
vxe=var(x,1);
vye=var(y,1);
Kreuzkovarianz (für streng periodische Signale keine Bessel-Korrektur, da Mittelwert keine Schätzgröße, ohne Rauschen exakt, mit Rauschen erwartungstreu) cyx=cxy*=1N⋅∑i=0N−1(xi−x‾)*(yi−y‾)=(1N⋅∑i=0N−1xi*yi)−x‾*y‾
cov(y,x,ddof=0)[0,1]
sum(conj(x).*y)/N-conj(mean(x))*mean(y)
(Kreuz-)Korrelationskoeffizient (für streng periodische Signale keine Bessel-Korrektur, da Mittelwert keine Schätzgröße, ohne Rauschen exakt, mit Rauschen systematische Fehler) γyx=γxy*=cyxsx2sy2=∑i=0N−1(xi−x‾)*(yi−y‾)(∑i=0N−1|xi−x‾|2)(∑i=0N−1|yi−y‾|2)
=(∑i=0N−1xi*yi)−Nx‾*y‾[(∑i=0N−1|xi|2)−N|x‾|2][(∑i=0N−1|yi|2)−N|y‾|2]
corrcoef(y,x,ddof=0)[0,1]
(sum(conj(x).*y)/N-conj(mean(x))*mean(y))/sqrt(var(x,1)*var(y,1))
Kreuzkorrelationsfunktion (ohne Rauschen exakt, mit Rauschen erwartungstreu) Ryx,k=Ryx(τk)=Rxy*(−τk)=1N⋅∑i=0N−1xi*y(i+k)modN
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
Ryx=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])*y[(i+k)%N]
  
  Ryx[k]=sum1/float(N)

plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
show()
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
tau=arange(0,N)*dt
X=fft(x)
Y=fft(y)
Pyx=conj(X)*Y/N**2
Ryx=N*ifft(Pyx)
plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
show()
 direkte Berechnung:
Ryx=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))*y(mod(i+k-1,N)+1);
end
Ryx(k+1)=sum1/N;
end
plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
tau=(0:N-1)*dt;
X=fft(x);
Y=fft(y);
Pyx=conj(X).*Y/N^2;
Ryx=N*ifft(Pyx);
plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
Berechnung für den Fall, dass die Signallängen Nx bzw. Ny von der anzuwendende Periode des Signals Ñ abweichen (Ñ≤Nx, Ñ≤Ny) Ryx,k=Ryx(τk)=Rxy*(−τk)=∑i=0Nx−1∑j=0Ny−1j−i≡kmodÑxi*yj∑i=0Nx−1∑j=0Ny−1j−i≡kmodÑ1
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
NT=60 #Periodenlänge NT≤N
Nx=len(x)
Ny=len(y)
R1=zeros(NT)+0j
R0=zeros(NT)
for i in range(0,Nx):
  for j in range(0,Ny):
    R1[(j-i)%NT]+=conj(x[i])*y[j]
    R0[(j-i)%NT]+=1
  

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

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

X1=fft(xp)
X0=fft(np)
Ny=len(x)
yp=zeros(NT)+0j
np=zeros(NT)
for i in range(0,Ny):
  yp[i%NT]+=y[i]
  np[i%NT]+=1

Y1=fft(yp)
Y0=fft(np)
tau=arange(0,NT)*dt
P1=conj(X1)*Y1/NT**2
P0=conj(X0)*Y0/NT**2
Ryx=ifft(P1)/real(ifft(P0))
plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
show()
 direkte Berechnung:
NT=60; %Periodenlänge NT≤N
Nx=length(x);
Ny=length(y);
R1=zeros(1,NT)+0j;
R0=zeros(1,NT);
for i=1:Nx
for j=1:Ny
R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+conj(x(i))*y(j);
R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+1;
end
end
tau=(0:NT-1)*dt;
Ryx=zeros(1,NT)+0j;
for k=1:NT
Ryx(k)=R1(k)/R0(k);
end
plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
NT=60; %Periodenlänge NT≤N
Nx=length(x);
Ny=length(y);
xp=zeros(1,NT)+0j;
np=zeros(1,NT);
for i=1:Nx
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);
yp=zeros(1,NT)+0j;
np=zeros(1,NT);
for i=1:Ny
yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+y(i);
np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1;
end
Y1=fft(yp);
Y0=fft(np);
tau=(0:NT-1)*dt;
P1=conj(X1).*Y1/NT^2;
P0=conj(X0).*Y0/NT^2;
Ryx=ifft(P1)./real(ifft(P0));
plot(tau,real(Ryx),'o',tau,imag(Ryx),'o')
Kreuzkovarianzfunktion (ohne Rauschen exakt, mit Rauschen erwartungstreu) Cyx,k=Cyx(τk)=Cxy*(−τk)=1N⋅∑i=0N−1(xi−x‾)*[y(i+k)modN−y‾]
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
Cyx=zeros(N)+0j
tau=arange(0,N)*dt
mxe=mean(x)
mye=mean(y)
for k in range(0,N):
  sum1=0+0j
  for i in range(0,N):
    sum1+=conj(x[i]-mxe)*(y[(i+k)%N]-mye)
  
  Cyx[k]=sum1/float(N)

plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
show()
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
tau=arange(0,N)*dt
mxe=mean(x)
mye=mean(y)
X=fft(x-mxe)
Y=fft(y-mye)
Pyx=conj(X)*Y/N**2
Cyx=N*ifft(Pyx)
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
show()
 direkte Berechnung:
Cyx=zeros(1,N)+0j;
tau=(0:N-1)*dt;
mxe=mean(x);
mye=mean(y);
for k=0:N-1
sum1=0+0j;
for i=1:N
sum1=sum1+conj(x(i)-mxe)*(y(mod(i+k-1,N)+1)-mye);
end
Cyx(k+1)=sum1/N;
end
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
tau=(0:N-1)*dt;
mxe=mean(x);
mye=mean(y);
X=fft(x-mxe);
Y=fft(y-mye);
Pyx=conj(X).*Y/N^2;
Cyx=N*ifft(Pyx);
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Berechnung für den Fall, dass die Signallängen Nx bzw. Ny von der anzuwendende Periode des Signals Ñ abweichen (Ñ≤Nx, Ñ≤Ny) Cyx,k=Cyx(τk)=Cxy*(−τk)=∑i=0Nx−1∑j=0Ny−1j−i≡kmodÑ(xi−x‾)*(yj−y‾)∑i=0Nx−1∑j=0Ny−1j−i≡kmodÑ1
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
mxe=mean(x)
mye=mean(y)
NT=60 #Periodenlänge NT≤N
Nx=len(x)
Ny=len(y)
C1=zeros(NT)+0j
C0=zeros(NT)
for i in range(0,Nx):
  for j in range(0,Ny):
    C1[(j-i)%NT]+=conj(x[i]-mxe)*(y[j]-mye)
    C0[(j-i)%NT]+=1
  

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

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

X1=fft(xp)
X0=fft(np)
yp=zeros(NT)+0j
np=zeros(NT)
for i in range(0,Ny):
  yp[i%NT]+=y[i]-mye
  np[i%NT]+=1

Y1=fft(yp)
Y0=fft(np)
tau=arange(0,NT)*dt
P1=conj(X1)*Y1/NT**2
P0=conj(X0)*Y0/NT**2
Cyx=ifft(P1)/real(ifft(P0))
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
show()
 direkte Berechnung:
mxe=mean(x);
mye=mean(y);
NT=60; %Periodenlänge NT≤N
Nx=length(x);
Ny=length(y);
C1=zeros(1,NT)+0j;
C0=zeros(1,NT);
for i=1:Nx
for j=1:Ny
C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+conj(x(i)-mxe)*(y(j)-mye);
C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+1;
end
end
tau=(0:NT-1)*dt;
Cyx=zeros(1,NT)+0j;
for k=1:NT
Cyx(k)=C1(k)/C0(k);
end
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
mxe=mean(x);
mye=mean(y);
NT=60; %Periodenlänge NT≤N
Nx=length(x);
Ny=length(y);
xp=zeros(1,NT)+0j;
np=zeros(1,NT);
for i=1:Nx
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
X1=fft(xp);
X0=fft(np);
yp=zeros(1,NT)+0j;
np=zeros(1,NT);
for i=1:Ny
yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+y(i)-mye;
np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1;
end
Y1=fft(yp);
Y0=fft(np);
tau=(0:NT-1)*dt;
P1=conj(X1).*Y1/NT^2;
P0=conj(X0).*Y0/NT^2;
Cyx=ifft(P1)./real(ifft(P0));
plot(tau,real(Cyx),'o',tau,imag(Cyx),'o')
Kreuzkorrelationskoeffizientenfunktion (ohne Rauschen exakt, mit Rauschen systematische Fehler) ρyx,k=ρyx(τk)=ρxy*(−τk)=Cyx,ksx2sy2=1Nsx2sy2⋅∑i=0N−1(xi−x‾)*[y(i+k)modN−y‾]
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
rhoyx=zeros(N)+0j
tau=arange(0,N)*dt
mxe=mean(x)
mye=mean(y)
vxe=var(x)
vye=var(y)
for k in range(0,N):
  sum1=0+0j
  for i in range(0,N):
    sum1+=conj(x[i]-mxe)*(y[(i+k)%N]-mye)
  
  rhoyx[k]=sum1/float(N)/sqrt(vxe*vye)

plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
show()
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
tau=arange(0,N)*dt
mxe=mean(x)
mye=mean(y)
vxe=var(x)
vye=var(y)
X=fft(x-mxe)
Y=fft(y-mye)
Pyx=conj(X)*Y/N**2
rhoyx=N*ifft(Pyx)/sqrt(vxe*vye)
plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
show()
 direkte Berechnung:
rhoyx=zeros(1,N)+0j;
tau=(0:N-1)*dt;
mxe=mean(x);
mye=mean(y);
vxe=var(x,1);
vye=var(y,1);
for k=0:N-1
sum1=0+0j;
for i=1:N
sum1=sum1+conj(x(i)-mxe)*(y(mod(i+k-1,N)+1)-mye);
end
rhoyx(k+1)=sum1/N/sqrt(vxe*vye);
end
plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
tau=(0:N-1)*dt;
mxe=mean(x);
mye=mean(y);
vxe=var(x,1);
vye=var(y,1);
X=fft(x-mxe);
Y=fft(y-mye);
Pyx=conj(X).*Y/N^2;
rhoyx=N*ifft(Pyx)/sqrt(vxe*vye);
plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
Berechnung für den Fall, dass die Signallängen Nx bzw. Ny von der anzuwendende Periode des Signals Ñ abweichen (Ñ≤Nx, Ñ≤Ny) ρyx,k=ρyx(τk)=ρxy*(−τk)=Cyx,ksx2sy2=∑i=0Nx−1∑j=0Ny−1j−i≡kmodÑ(xi−x‾)*(yj−y‾)sx2sy2⋅∑i=0Nx−1∑j=0Ny−1j−i≡kmodÑ1
τk=kΔt
k=−∞…∞ 		direkte Berechnung:
mxe=mean(x)
mye=mean(y)
vxe=var(x)
vye=var(y)
NT=60 #Periodenlänge NT≤N
Nx=len(x)
Ny=len(y)
C1=zeros(NT)+0j
C0=zeros(NT)
for i in range(0,Nx):
  for j in range(0,Ny):
    C1[(j-i)%NT]+=conj(x[i]-mxe)*(y[j]-mye)
    C0[(j-i)%NT]+=1
  

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

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

X1=fft(xp)
X0=fft(np)
yp=zeros(NT)+0j
np=zeros(NT)
for i in range(0,Ny):
  yp[i%NT]+=y[i]-mye
  np[i%NT]+=1

Y1=fft(yp)
Y0=fft(np)
tau=arange(0,NT)*dt
P1=conj(X1)*Y1/NT**2
P0=conj(X0)*Y0/NT**2
rhoyx=ifft(P1)/real(ifft(P0))/sqrt(vxe*vye)
plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
show()
 direkte Berechnung:
mxe=mean(x);
mye=mean(y);
vxe=var(x,1);
vye=var(y,1);
NT=60; %Periodenlänge NT≤N
Nx=length(x);
Ny=length(y);
C1=zeros(1,NT)+0j;
C0=zeros(1,NT);
for i=1:Nx
for j=1:Ny
C1(mod(j-i,NT)+1)=C1(mod(j-i,NT)+1)+conj(x(i)-mxe)*(y(j)-mye);
C0(mod(j-i,NT)+1)=C0(mod(j-i,NT)+1)+1;
end
end
tau=(0:NT-1)*dt;
rhoyx=zeros(1,NT)+0j;
for k=1:NT
rhoyx(k)=C1(k)/C0(k)/sqrt(vxe*vye);
end
plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
Berechnung aus Leistungsspektrum mittels Wiener-Chintschin-Theorem:
mxe=mean(x);
mye=mean(y);
vxe=var(x,1);
vye=var(y,1);
NT=60; %Periodenlänge NT≤N
Nx=length(x);
Ny=length(y);
xp=zeros(1,NT)+0j;
np=zeros(1,NT);
for i=1:Nx
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
X1=fft(xp);
X0=fft(np);
yp=zeros(1,NT)+0j;
np=zeros(1,NT);
for i=1:Ny
yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+y(i)-mye;
np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1;
end
Y1=fft(yp);
Y0=fft(np);
tau=(0:NT-1)*dt;
P1=conj(X1).*Y1/NT^2;
P0=conj(X0).*Y0/NT^2;
rhoyx=ifft(P1)./real(ifft(P0))/sqrt(vxe*vye);
plot(tau,real(rhoyx),'o',tau,imag(rhoyx),'o')
Kreuzleistungsspektrum (ohne Rauschen exakt, mit Rauschen erwartungstreu) Xj=FFT{xi}=∑i=0N−1xi𝐞−2π𝐢ij/N
Yj=FFT{yi}=∑i=0N−1yi𝐞−2π𝐢ij/N
(imaginäre Einheit 𝐢)
Pyx,j=Pyx(fj)=Xj*YjN2
fj=jNΔt
j=−⌊N/2⌋…⌊(N−1)/2⌋ 		direkte Berechnung:
from numpy.fft import *
f=roll(arange(-(N//2),N-(N//2))/float(N*dt),N-(N//2))
X=fft(x)
Y=fft(y)
Pyx=conj(X)*Y/N**2
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
show()
Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
from numpy.fft import *
Ryx=zeros(N)+0j
for k in range(0,N):
  sum1=0+0j
  for i in range(0,N):
    sum1+=conj(x[i])*y[(i+k)%N]
  
  Ryx[k]=sum1/N

f=roll(arange(-(N//2),N-(N//2))/float(N*dt),N-(N//2))
Pyx=fft(Ryx)/float(N)
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
show()
 direkte Berechnung:
f=circshift(((-N/2:N-N/2-1))/(N*dt),[0;N/2]);
X=fft(x);
Y=fft(y);
Pyx=conj(X).*Y/N^2;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Berechnung aus Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
Ryx=zeros(1,N)+0j;
for k=0:N-1
sum1=0+0j;
for i=1:N
sum1=sum1+conj(x(i))*y(mod(i+k-1,N)+1);
end
Ryx(k+1)=sum1/N;
end
f=circshift(((-N/2:N-N/2-1))/(N*dt),[0;N/2]);
Pyx=fft(Ryx)/N;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Berechnung für den Fall, dass die Signallängen Nx bzw. Ny von der anzuwendende Periode des Signals Ñ abweichen (Ñ≤Nx, Ñ≤Ny) (wegen der nötigen Entfaltung mit dem Leistungsspektrum der Überlappungen nur über die Korrelationsfunktion)
Pyx,j=Pyx(fj)=1Ñ⋅FFT{Ryx,k}=1Ñ⋅∑k=0Ñ−1Ryx,k𝐞−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
Nx=len(x)
Ny=len(y)
xp=zeros(NT)+0j
np=zeros(NT)
for i in range(0,Nx):
  xp[i%NT]+=x[i]
  np[i%NT]+=1

X1=fft(xp)
X0=fft(np)
yp=zeros(NT)+0j
np=zeros(NT)
for i in range(0,Ny):
  yp[i%NT]+=y[i]
  np[i%NT]+=1

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

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

f=roll(arange(-(NT//2),(NT+1)//2)/float(NT*dt),(NT+1)//2)
Pyx=fft(Ryx)/float(NT)
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
show()
 Berechnung aus Korrelationsfunktion mittels zweifacher Anwendung des Wiener-Chintschin-Theorems, direkte Bestimmung der primären Spektren:
NT=60; %Periodenlänge NT≤N
Nx=length(x);
Ny=length(y);
xp=zeros(1,NT)+0j;
np=zeros(1,NT);
for i=1:Nx
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);
yp=zeros(1,NT)+0j;
np=zeros(1,NT);
for i=1:Ny
yp(mod(i-1,NT)+1)=yp(mod(i-1,NT)+1)+y(i);
np(mod(i-1,NT)+1)=np(mod(i-1,NT)+1)+1;
end
Y1=fft(yp);
Y0=fft(np);
P1=conj(X1).*Y1/NT^2;
P0=conj(X0).*Y0/NT^2;
Ryx=ifft(P1)./real(ifft(P0));
f=circshift((-fix(NT/2):fix((NT-1)/2))/(NT*dt),[0;fix((NT+1)/2)]);
Pyx=fft(Ryx)/NT;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')
Berechnung aus direkt bestimmter Korrelationsfunktion mittels Wiener-Chintschin-Theorem:
NT=60; %Periodenlänge NT≤N
Nx=length(x);
Ny=length(y);
R1=zeros(1,NT)+0j;
R0=zeros(1,NT);
for i=1:Nx
for j=1:Ny
R1(mod(j-i,NT)+1)=R1(mod(j-i,NT)+1)+conj(x(i))*y(j);
R0(mod(j-i,NT)+1)=R0(mod(j-i,NT)+1)+1;
end
end
Ryx=zeros(1,NT)+0j;
for k=1:NT
Ryx(k)=R1(k)/R0(k);
end
f=circshift((-fix(NT/2):fix((NT-1)/2))/(NT*dt),[0;fix((NT+1)/2)]);
Pyx=fft(Ryx)/NT;
plot(f,real(Pyx),'o',f,imag(Pyx),'o')