Signal- und Messdatenverarbeitung

Interpolation, Extrapolation und Regression

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
lineare Interpolation und -extrapolation (x,y): Vektoren mit zwei zu interpolierenden Wertepaaren
c=(c1c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte der linearen Interpolation.
c=polyfit(x,y,1)
def f(x):
  return polyval(c,x)
c=polyfit(x,y,1);
f=@(x)(polyval(c,x));
quadratische Interpolation und -extrapolation (mit Polynom zweiten Grades) (x,y): Vektoren mit drei zu interpolierenden Wertepaaren
c=(c2c1c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte der linearen Interpolation.
c=polyfit(x,y,2)
def f(x):
  return polyval(c,x)
c=polyfit(x,y,2);
f=@(x)(polyval(c,x));
Polynominterpolation und -extrapolation (x,y): Vektoren mit J zu interpolierenden Wertepaaren
c=(cJ−1cJ−2⋯c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte des Interpolationspolynoms.
c=polyfit(x,y,len(x)-1)
def f(x):
  return polyval(c,x)
c=polyfit(x,y,length(x)-1);
f=@(x)(polyval(c,x));
Polynomregression (ohne Gewichtung) (x,y): Vektoren mit J Wertepaaren
K: Grad des anzupassenden Polynoms mit K<J−1
c=(cKcK−1⋯c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte des angepassten Polynoms.
c=polyfit(x,y,K)
def f(x):
  return polyval(c,x)
c=polyfit(x,y,K);
f=@(x)(polyval(c,x));
Polynomregression mit Gewichtung (x,y): Vektoren mit J Wertepaaren
w: Vektor mit J Gewichten
K: Grad des anzupassenden Polynoms mit K<J−1 c=(cKcK−1⋯c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte des angepassten Polynoms.
from numpy.linalg import *
cm=zeros([K+1,K+1])
cs=zeros(K+1)
for i in range(0,K+1):
  for j in range(0,K+1):
    cm[i,j]=sum(array(w)*array(x)**(i+j))
  cs[i]=sum(array(w)*array(y)*array(x)**i)

cc=solve(cm,cs)
def f(x):
  return polyval(flipud(cc),x)
cm=zeros(K+1,K+1);
cs=zeros(K+1,1);
for i=1:K+1
for j=1:K+1
cm(i,j)=sum(w.*x.^(i+j-2));
end
cs(i)=sum(w.*y.*x.^(i-1));
end
cc=linsolve(cm,cs);
f=@(x)(polyval(flipud(cc),x));
Polynomregression (ohne Gewichtung) unter der Randbedingung eines verschwindenden Absolutglieds (x,y): Vektoren mit J Wertepaaren
K: Grad des anzupassenden Polynoms mit K<J−1 c=(cKcK−1⋯c0=0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte des angepassten Polynoms.
from numpy.linalg import *
cm=zeros([K,K])
cs=zeros(K)
for i in range(0,K):
  for j in range(0,K):
    cm[i,j]=sum(array(x)**(i+j+2))
  cs[i]=sum(array(y)*array(x)**(i+1))

cc=concatenate(([0],solve(cm,cs)))
def f(x):
  return polyval(flipud(cc),x)
cm=zeros(K,K);
cs=zeros(K,1);
for i=1:K
for j=1:K
cm(i,j)=sum(x.^(i+j));
end
cs(i)=sum(y.*x.^i);
end
cc=[0;linsolve(cm,cs)];
f=@(x)(polyval(flipud(cc),x));
Polynomregression mit Gewichtung unter der Randbedingung eines verschwindenden Absolutglieds (x,y): Vektoren mit J Wertepaaren
w: Vektor mit J Gewichten
K: Grad des anzupassenden Polynoms mit K<J−1 c=(cKcK−1⋯c0=0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte des angepassten Polynoms.
from numpy.linalg import *
cm=zeros([K,K])
cs=zeros(K)
for i in range(0,K):
  for j in range(0,K):
    cm[i,j]=sum(array(w)*array(x)**(i+j+2))
  cs[i]=sum(array(w)*array(y)*array(x)**(i+1))

cc=concatenate(([0],solve(cm,cs)))
def f(x):
  return polyval(flipud(cc),x)
cm=zeros(K,K);
cs=zeros(K,1);
for i=1:K
for j=1:K
cm(i,j)=sum(w.*x.^(i+j));
end
cs(i)=sum(w.*y.*x.^i);
end
cc=[0;linsolve(cm,cs)];
f=@(x)(polyval(flipud(cc),x));
Polynominterpolation und -extrapolation unter der Randbedingung der Symmetrie zur y-Achse (x,y): Vektoren mit J zu interpolierenden Wertepaaren
c=(c2J−2c2J−4⋯c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte des Interpolationspolynoms.
c=polyfit(array(x)**2,y,len(x)-1)
def f(x):
  return polyval(c,array(x)**2)
c=polyfit(x.^2,y,length(x)-1);
f=@(x)(polyval(c,x.^2));
Polynomregression (ohne Gewichtung) unter der Randbedingung der Symmetrie zur y-Achse (x,y): Vektoren mit J Wertepaaren
K: Grad des anzupassenden Polynoms (gerade) mit K<2J−2
c=(cKcK−2⋯c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte des angepassten Polynoms.
from numpy.linalg import *
cm=zeros([K//2+1,K//2+1])
cs=zeros(K//2+1)
for i in range(0,K//2+1):
  for j in range(0,K//2+1):
    cm[i,j]=sum(array(x)**(2*(i+j)))
  cs[i]=sum(array(y)*array(x)**(2*i))

cc=solve(cm,cs)
def f(x):
  return polyval(flipud(cc),x**2)
cm=zeros(fix(K/2)+1,fix(K/2)+1);
cs=zeros(fix(K/2)+1,1);
for i=1:fix(K/2)+1
for j=1:fix(K/2)+1
cm(i,j)=sum(x.^(2*(i+j-2)));
end
cs(i)=sum(y.*x.^(2*i-2));
end
cc=linsolve(cm,cs);
f=@(x)(polyval(flipud(cc),x.^2));
Polynomregression mit Gewichtung unter der Randbedingung der Symmetrie zur y-Achse (x,y): Vektoren mit J Wertepaaren
w: Vektor mit J Gewichten
K: Grad des anzupassenden Polynoms (gerade) mit K<2J−2
c=(cKcK−2⋯c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte des angepassten Polynoms.
from numpy.linalg import *
cm=zeros([K//2+1,K//2+1])
cs=zeros(K//2+1)
for i in range(0,K//2+1):
  for j in range(0,K//2+1):
    cm[i,j]=sum(array(x)**(2*(i+j)))
  cs[i]=sum(array(y)*array(x)**(2*i))

cc=solve(cm,cs)
def f(x):
  return polyval(flipud(cc),x**2)
cm=zeros(fix(K/2)+1,fix(K/2)+1);
cs=zeros(fix(K/2)+1,1);
for i=1:fix(K/2)+1
for j=1:fix(K/2)+1
cm(i,j)=sum(w.*x.^(2*(i+j-2)));
end
cs(i)=sum(w.*y.*x.^(2*i-2));
end
cc=linsolve(cm,cs);
f=@(x)(polyval(flipud(cc),x.^2));
exponentielle Interpolation und -extrapolation (x,y): Vektoren mit zwei zu interpolierenden Wertepaaren
c=(c1c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte der exponentiellen Interpolation.
c=polyfit(-array(x),log(y),1)
def f(x):
  return exp(c[1]-c[0]*x)
c=polyfit(-x,log(y),1);
f=@(x)(exp(c(2)-c(1)*x));
exponentielle Regression ohne Gewichtung (x,y): Vektoren mit J>2 Wertepaaren
c=(c1c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte der exponentiellen Interpolation.
c=polyfit(-array(x),log(y),1)
def f(x):
  return exp(c[1]-c[0]*x)
c=polyfit(-x,log(y),1);
f=@(x)(exp(c(2)-c(1)*x));
exponentielle Regression mit Gewichtung (x,y): Vektoren mit J>2 Wertepaaren
w: Vektor mit J Gewichten
c=(c1c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte der exponentiellen Interpolation.
from numpy.linalg import *
cm=zeros([2,2])
cs=zeros(2)
for i in range(0,2):
  for j in range(0,2):
    cm[i,j]=sum(array(w)*array(x)**(i+j))
  cs[i]=sum(array(w)*array(log(y))*array(x)**i)

cc=solve(cm,cs)
def f(x):
  return exp(cc[0])*exp(cc[1]*x)
cm=zeros(2,2);
cs=zeros(2,1);
for i=1:2
for j=1:2
cm(i,j)=sum(w.*x.^(i+j-2));
end
cs(i)=sum(w.*log(y).*x.^(i-1));
end
cc=linsolve(cm,cs);
f=@(x)(exp(cc(1))*exp(cc(2)*x));
Gauß-Interpolation und -extrapolation (x,y): Vektoren mit drei zu interpolierenden Wertepaaren
c=(c2c1c0): Koeffizienten-(zeilen-)vektor
Die Funktion f(x) liefert die Funktionswerte der linearen Interpolation.
c=polyfit(x,log(array(y)),2)
def f(x):
  return exp(polyval(c,x))
c=polyfit(x,log(y),2);
f=@(x)(exp(polyval(c,x)));
Gauß-Interpolation der drei Werte um das Maximum in einem Vektor und Bestimmung der Position des Maximums in der interpolierten Funktion y: Vektor mit Werten
X: Argument des Maximums (Index)
δ: Verschiebung des Maximums der interpolierten Gauß-Funktion gegenüber dem Index des maximalen Wertes
X=maximum(1,minimum(len(y)-2,argmax(y)))
delta=(log(y[X-1])-log(y[X+1]))/(2*(log(y[X+1])-2*log(y[X])+log(y[X-1])))
[Y,X]=max(y);
X=max(2,min(length(y)-1,X));
delta=(log(y(X-1))-log(y(X+1)))/(2*(log(y(X+1))-2*log(y(X))+log(y(X-1))));
(stückweise) Sample-and-Hold-Interpolation (xp,yp): Eingangswertepaare in aufsteigender Reihenfolge ihrer Argumente
xq: Vektor mit Argumenten der interpolierten Funktion
yq: Ergebnisvektor mit Werten der interpolierten Funktion
from scipy.interpolate import *
f=interp1d(xp,yp,'zero')
yq=f(array(xq))
yq=zeros(size(xq));
for q=1:length(xq)
yq(q)=yp(find(xp<=xq(q),1,'last'));
end
(stückweise) Nearest-Neighbor-Interpolation (xp,yp): Eingangswertepaare in aufsteigender Reihenfolge ihrer Argumente
xq: Vektor mit Argumenten der interpolierten Funktion
yq: Ergebnisvektor mit Werten der interpolierten Funktion
from scipy.interpolate import *
f=interp1d(xp,yp,'nearest')
yq=f(array(xq))
yq=interp1(xp,yp,xq,'nearest');
stückweise lineare Interpolation (xp,yp): Eingangswertepaare in aufsteigender Reihenfolge ihrer Argumente
xq: Vektor mit Argumenten der interpolierten Funktion
yq: Ergebnisvektor mit Werten der interpolierten Funktion
from scipy.interpolate import *
f=interp1d(xp,yp,'linear')
yq=f(array(xq))

oder

from scipy.interpolate import *
f=interp1d(xp,yp)
yq=f(array(xq))
yq=interp1(xp,yp,xq,'linear');

oder

yq=interp1(xp,yp,xq);
(stückweise) Interpolation mit natürlichen Splines (xp,yp): Eingangswertepaare in aufsteigender Reihenfolge ihrer Argumente
xq: Vektor mit Argumenten der interpolierten Funktion
yq: Ergebnisvektor mit Werten der interpolierten Funktion
from numpy.linalg import *
cm=zeros([len(xp)-2,len(xp)-2])
for i in range(0,len(xp)-2):
  cm[i,i]=(xp[i+2]-xp[i])/3.0

for i in range(0,len(xp)-3):
  cm[i,i+1]=(xp[i+2]-xp[i+1])/6.0
  cm[i+1,i]=(xp[i+2]-xp[i+1])/6.0

cs=zeros(len(xp)-2)
for i in range(0,len(xp)-2):
  cs[i]=float(yp[i+2]-yp[i+1])/(xp[i+2]-xp[i+1])-float(yp[i+1]-yp[i])/(xp[i+1]-xp[i])

cc=solve(cm,cs)
yq=zeros(size(xq))
nextLowest=lambda seq,x: min([(x-i,i) for i in seq if x>=i] or [(0,None)])[1]
for q in range(0,len(xq)):
  i=minimum(len(xp)-2,xp.index(nextLowest(xp,xq[q])))
  yq[q]=float(xp[i+1]-xq[q])/(xp[i+1]-xp[i])*yp[i]+float(xq[q]-xp[i])/(xp[i+1]-xp[i])*yp[i+1]
  if i>0:
    yq[q]+=cc[i-1]*(xp[i+1]-xp[i])*(xp[i+1]-xq[q])/6.0*((float(xp[i+1]-xq[q])/(xp[i+1]-xp[i]))**2-1)
  if i<len(xp)-2:
    yq[q]+=cc[i]*(xp[i+1]-xp[i])*(xq[q]-xp[i])/6.0*((float(xq[q]-xp[i])/(xp[i+1]-xp[i]))**2-1)
  
cm=zeros(length(xp)-2,length(xp)-2);
for i=1:length(xp)-2
cm(i,i)=(xp(i+2)-xp(i))/3;
end
for i=1:length(xp)-3
cm(i,i+1)=(xp(i+2)-xp(i+1))/6;
cm(i+1,i)=(xp(i+2)-xp(i+1))/6;
end
cs=zeros(length(xp)-2,1);
for i=1:length(xp)-2
cs(i)=(yp(i+2)-yp(i+1))/(xp(i+2)-xp(i+1))-(yp(i+1)-yp(i))/(xp(i+1)-xp(i));
end
cc=linsolve(cm,cs);
yq=zeros(size(xq));
for q=1:length(xq)
i=min(length(xp)-1,find(xp<=xq(q),1,'last'));
yq(q)=(xp(i+1)-xq(q))/(xp(i+1)-xp(i))*yp(i)+(xq(q)-xp(i))/(xp(i+1)-xp(i))*yp(i+1);
if i>1
yq(q)=yq(q)+cc(i-1)*(xp(i+1)-xp(i))*(xp(i+1)-xq(q))/6*(((xp(i+1)-xq(q))/(xp(i+1)-xp(i)))^2-1);
end
if i<length(xp)-1
yq(q)=yq(q)+cc(i)*(xp(i+1)-xp(i))*(xq(q)-xp(i))/6*(((xq(q)-xp(i))/(xp(i+1)-xp(i)))^2-1);
end
end
(stückweise) Interpolation mit natürlichen Splines mit äquidistanten Abtaststellen yp: Vektor mit n Eingangswerten
xq: Vektor mit Argumenten der interpolierten Funktion (Python: 0…n−1, Matlab: 1…n)
yq: Ergebnisvektor mit Werten der interpolierten Funktion
cm=zeros([len(yp)-2,len(yp)-2])
for i in range(0,len(yp)-2):
  cm[i,i]=2.0/3

for i in range(0,len(yp)-3):
  cm[i,i+1]=1.0/6
  cm[i+1,i]=1.0/6

cs=zeros(len(yp)-2)
for i in range(0,len(yp)-2):
  cs[i]=yp[i]-2*yp[i+1]+yp[i+2]

cc=solve(cm,cs)
yq=zeros(size(xq))
for q in range(0,len(xq)):
  i=minimum(len(yp)-2,int(floor(xq[q])))
  yq[q]=(i+1-xq[q])*yp[i]+(xq[q]-i)*yp[i+1]
  if i>0:
    yq[q]+=cc[i-1]*(i+1-xq[q])*((i+1-xq[q])**2-1)/6.0
  if i<len(yp)-2:
    yq[q]+=cc[i]*(xq[q]-i)*((xq[q]-i)**2-1)/6.0
  
cm=zeros(length(yp)-2,length(yp)-2);
for i=1:length(yp)-2
cm(i,i)=2/3;
end
for i=1:length(yp)-3
cm(i,i+1)=1/6;
cm(i+1,i)=1/6;
end
cs=zeros(length(yp)-2,1);
for i=1:length(yp)-2
cs(i)=yp(i)-2*yp(i+1)+yp(i+2);
end
cc=linsolve(cm,cs);
yq=zeros(size(xq));
for q=1:length(xq)
i=min(length(yp)-1,floor(xq(q)));
yq(q)=(i+1-xq(q))*yp(i)+(xq(q)-i)*yp(i+1);
if i>1
yq(q)=yq(q)+cc(i-1)*(i+1-xq(q))/6*((i+1-xq(q))^2-1);
end
if i<length(yp)-1
yq(q)=yq(q)+cc(i)*(xq(q)-i)/6*((xq(q)-i)^2-1);
end
end
(stückweise) Splineinterpolation mit not-a-knot-Bedingungen (xp,yp): Eingangswertepaare in aufsteigender Reihenfolge ihrer Argumente
xq: Vektor mit Argumenten der interpolierten Funktion
yq: Ergebnisvektor mit Werten der interpolierten Funktion
from scipy.interpolate import *
yq=splev(xq,splrep(xp,yp,s=0.0))
#Für die SciPy-Standardprozedur interp1d(xp,yp,'cubic') sind die verwendeten Randbedingungen unklar.
yq=spline(xp,yp,xq);
(stückweise) Splineinterpolation mit periodischen Randbedingungen (xp,yp): Eingangswertepaare in aufsteigender Reihenfolge ihrer Argumente
xq: Vektor mit Argumenten der interpolierten Funktion
yq: Ergebnisvektor mit Werten der interpolierten Funktion
from scipy.interpolate import *
yq=splev(xq,splrep(xp,yp,s=0.0,per=1))
 fehlt in Matlab/Octave
parametrische, stückweise lineare Interpolation (xp,yp): Vektoren mit zu interpolierenden Wertepaaren
r: Vektor mit Argumenten der Parameterfunktion (Python: 0…n−1, Matlab: 1…n)
(xr,yr): Ergebnisvektoren mit interpolierten Wertepaaren
from scipy.interpolate import *
fx=interp1d(arange(0,len(xp)),xp)
fy=interp1d(arange(0,len(yp)),yp)
xr=fx(r)
yr=fy(r)
xr=interp1(xp,r);
yr=interp1(yp,r);
parametrische Interpolation mit (stückweisen) Splines mit not-a-knot-Bedingungen (xp,yp): Vektoren mit zu interpolierenden Wertepaaren
r: Vektor mit Argumenten der Parameterfunktion (Python: 0…n−1, Matlab: 1…n)
(xr,yr): Ergebnisvektoren mit interpolierten Wertepaaren
from scipy.interpolate import *
xr=splev(r,splrep(arange(0,len(xp)),xp,s=0.0))
yr=splev(r,splrep(arange(0,len(yp)),yp,s=0.0))
xr=spline((1:length(xp)),xp,r);
yr=spline((1:length(yp)),yp,r);
Whittaker-Shannon-Interpolation (sinc-, si-Interpolation) mit äquidistanten Abtaststellen yp: Vektor mit n Eingangswerten an den Stellen 0…n−1 (Python) bzw. 1…n (Matlab)
xq: Vektor mit Argumenten der interpolierten Funktion (Python: 0…n−1, Matlab: 1…n)
yq: Ergebnisvektor mit Werten der interpolierten Funktion
yq=mean(yp)*ones(len(xq))
for i in range(0,size(yp)):
  yq+=(yp[i]-mean(yp))*sinc(xq-i)
yq=mean(yp)*ones(size(xq));
for i=1:length(yp)
yq=yq+(yp(i)-mean(yp)).*sinc(xq-i);
end
Whittaker-Shannon-Interpolation (sinc-, si-Interpolation) mit äquidistanten Abtaststellen für periodisch fortzusetzende Signale yp: Vektor mit n Eingangswerten an den Stellen 0…n−1 (Python) bzw. 1…n (Matlab)
xq: Vektor mit Argumenten der interpolierten Funktion (Python: 0…n−1, Matlab: 1…n)
yq: Ergebnisvektor mit Werten der interpolierten Funktion
def sinn(x,n):
  if sin(x)!=0:
    return sin(n*x)/sin(x)
  else:
    return n*cos(n*x)*cos(x)
  

yq=mean(yp)*ones(len(xq))
for i in range(0,len(xq)):
  if len(yp)%2==0:
    for j in range(0,len(yp)):
      yq[i]+=(yp[j]-mean(yp))*sinn(pi*(xq[i]-j)/float(len(yp)),len(yp))*cos(pi*(xq[i]-j)/float(len(yp)))/float(len(yp))
    
  else:
    for j in range(0,len(yp)):
      yq[i]+=(yp[j]-mean(yp))*sinn(pi*(xq[i]-j)/float(len(yp)),len(yp))/float(len(yp))
    
  
fehlt noch
zweidimensionale lineare Interpolation und -extrapolation (xp,yp,zp): Vektoren mit drei zu interpolierenden Wertetripeln
(xq,yq): Vektoren mit Argumenten der interpolierten Funktion
zq: Ergebnisvektor mit Werten der interpolierten Funktion
from numpy.linalg import *
c=solve(hstack((ones([3,1]),matrix(xp).T,matrix(yp).T)),matrix(zp).T)
zq=float(c[0])+float(c[1])*xq+float(c[2])*yq
c=linsolve([ones(length(xp),1) xp(1:3)' yp(1:3)'],zp(1:3)');
zq=c(1)+c(2)*xq+c(3)*yq;
zweidimensionale quadratische Interpolation und -extrapolation (xp,yp,zp): Vektoren mit sechs zu interpolierenden Wertetripeln
(xq,yq): Vektoren mit Argumenten der interpolierten Funktion
zq: Ergebnisvektor mit Werten der interpolierten Funktion
from numpy.linalg import *
c=solve(hstack((ones([6,1]),matrix(xp).T,matrix(yp).T,matrix(array(xp)**2).T,matrix(array(xp)*array(yp)).T,matrix(array(yp)**2).T)),matrix(zp).T)
zq=float(c[0])+float(c[1])*xq+float(c[2])*yq+float(c[3])*xq**2+float(c[4])*xq*yq+float(c[5])*yq**2
c=linsolve([ones(length(xp),1) xp(1:6)' yp(1:6)' xp(1:6)'.^2 xp(1:6)'.*yp(1:6)' yp(1:6)'.^2],zp(1:6)');
zq=c(1)+c(2)*xq+c(3)*yq+c(4)*xq.^2+c(5)*xq.*yq+c(6)*yq.^2;
stückweise bilineare Interpolation (eine nichtlineare Interpolation!) zp: Matrix mit zu interpolierenden Werten
(xq,yq): Vektoren mit Argumenten der interpolierten Funktion (Python: 0…n−1, Matlab: 1…n)
zq: Ergebnisvektor mit Werten der interpolierten Funktion
from scipy.interpolate import *
f=interp2d(arange(0,size(zp[0,:])),arange(0,size(zp[:,0])),zp,'linear')
zq=zeros(len(xq))
for i in range(0,len(xq)):
  zq[i]=f(xq[i],yq[i])

oder

from scipy.interpolate import *
f=interp2d(arange(0,size(zp[0,:])),arange(0,size(zp[:,0])),zp)
zq=zeros(len(xq))
for i in range(0,len(xq)):
  zq[i]=f(xq[i],yq[i])
zq=interp2(zp,xq,yq,'linear');

oder

zq=interp2(zp,xq,yq);
stückweise Interpolation mit bikubischen Splines mit not-a-knot-Bedingungen zp: Matrix mit zu interpolierenden Werten
(xq,yq): Vektoren mit Argumenten der interpolierten Funktion (Python: 0…n−1, Matlab: 1…n)
zq: Ergebnisvektor mit Werten der interpolierten Funktion
zq=zeros(len(xq))
for i in range(0,len(xq)):
  zs=zeros(len(zp[:,0]))
  for j in range(0,len(zp[:,0])):
    zs[j]=splev(xq[i],splrep(arange(0,len(zp[j,:])),zp[j,:],s=0.0))
  zq[i]=splev(yq[i],splrep(arange(0,len(zs)),zs,s=0.0))

#Weder für die SciPy-Standardprozedur interp2d(xp,yp,zp,'cubic') noch für bisplev(xq,yq,bisplrep(xp,yp,zp,s=0.0)) sind die verwendeten Randbedingungen klar.
#Die Interpolationsergebnisse dieser beiden Prozeduren sind darüber hinaus nicht separierbar.
zq=interp2(zp,xq,yq,'spline');
zweidimensionale Whittaker-Shannon-Interpolation (sinc-, si-Interpolation) mit äquidistanten Abtaststellen zp: Matrix mit zu interpolierenden Werten
(xq,yq): Vektoren mit Argumenten der interpolierten Funktion (Python: 0…n−1, Matlab: 1…n)
zq: Ergebnisvektor mit Werten der interpolierten Funktion
zq=mean(zp)*ones(len(xq))
for i in range(0,len(zp)):
  for j in range(0,len(zp[0])):
    zq+=(zp[i,j]-mean(zp))*sinc(yq-i)*sinc(xq-j)
  
zq=mean(mean(zp))*ones(size(xq));
[iz,jz]=size(zp);
for i=1:iz
for j=1:jz
zq=zq+(zp(i,j)-mean(mean(zp))).*sinc(yq-i).*sinc(xq-j);
end
end
zweidimensionale Whittaker-Shannon-Interpolation (sinc-, si-Interpolation) mit äquidistanten Abtaststellen für periodisch fortzusetzende Signale zp: Matrix mit zu interpolierenden Werten
(xq,yq): Vektoren mit Argumenten der interpolierten Funktion (Python: 0…n−1, Matlab: 1…n)
zq: Ergebnisvektor mit Werten der interpolierten Funktion
def sinn(x,n):
  if sin(x)!=0:
    return sin(n*x)/sin(x)
  else:
    return n*cos(n*x)*cos(x)
  

zq=mean(zp)*ones(len(xq))
for k in range(0,len(xq)):
  if (len(zp)%2==0) and (len(zp[0])%2==0):
    for i in range(0,len(zp)):
      for j in range(0,len(zp[0])):
        zq[k]+=(zp[i,j]-mean(zp))*sinn(pi*(xq[k]-j)/float(len(zp[0])),len(zp[0]))*cos(pi*(xq[k]-j)/float(len(zp[0])))/float(len(zp[0]))*sinn(pi*(yq[k]-i)/float(len(zp)),len(zp))*cos(pi*(yq[k]-i)/float(len(zp)))/float(len(zp))
      
    
  if (len(zp)%2==0) and (len(zp[0])%2==1):
    for i in range(0,len(zp)):
      for j in range(0,len(zp[0])):
        zq[k]+=(zp[i,j]-mean(zp))*sinn(pi*(xq[k]-j)/float(len(zp[0])),len(zp[0]))/float(len(zp[0]))*sinn(pi*(yq[k]-i)/float(len(zp)),len(zp))*cos(pi*(yq[k]-i)/float(len(zp)))/float(len(zp))
      
    
  if (len(zp)%2==1) and (len(zp[0])%2==0):
    for i in range(0,len(zp)):
      for j in range(0,len(zp[0])):
        zq[k]+=(zp[i,j]-mean(zp))*sinn(pi*(xq[k]-j)/float(len(zp[0])),len(zp[0]))*cos(pi*(xq[k]-j)/float(len(zp[0])))/float(len(zp[0]))*sinn(pi*(yq[k]-i)/float(len(zp)),len(zp))/float(len(zp))
      
    
  if (len(zp)%2==1) and (len(zp[0])%2==1):
    for i in range(0,len(zp)):
      for j in range(0,len(zp[0])):
        zq[k]+=(zp[i,j]-mean(zp))*sinn(pi*(xq[k]-j)/float(len(zp[0])),len(zp[0]))/float(len(zp[0]))*sinn(pi*(yq[k]-i)/float(len(zp)),len(zp))/float(len(zp))
      
    
  
fehlt noch