Signal- und Messdatenverarbeitung

Quellenanalyse

Empirical Mode Decomposition (EMD)

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
from numpy.linalg import *
Generieren eines Beispielsignals
N=1000
t=1.0*arange(0,N)
x=sin(0.05*t*(1+t/float(N)))+sin(0.02*t)+2*t**2/float(N**2)
plot(t,x)
show()
N=1000;
t=1*(0:N-1);
x=sin(0.05*t.*(1+t/N))+sin(0.02*t)+2*t.^2/N^2;
plot(t,x)
Definition der Spline-Interpolation (ähnlich den natürlichen Splines, aber mit einer an den Rändern des Interpolations- bzw. Extrapolationsbereiches verschwindenden zweiten Ableitung)
def mysplines(xp,yp,xleft,xright,xq):
  cm=zeros([4*(len(xp)-1),4*(len(xp)-1)])
  for i in range(0,len(xp)-1):
    cm[4*i,4*i]=0.0
    cm[4*i,4*i+1]=0.0
    cm[4*i,4*i+2]=1.0
    if i==0:
      cm[4*i,4*i+3]=3*xleft
    else:
      cm[4*i,4*i+3]=3*xp[i]
      cm[4*i,4*i-4]=0.0
      cm[4*i,4*i-3]=0.0
      cm[4*i,4*i-2]=-1.0
      cm[4*i,4*i-1]=-3*xp[i]
    cm[4*i+1,4*i]=1.0
    cm[4*i+1,4*i+1]=xp[i]
    cm[4*i+1,4*i+2]=xp[i]**2
    cm[4*i+1,4*i+3]=xp[i]**3
    cm[4*i+2,4*i]=1.0
    cm[4*i+2,4*i+1]=xp[i+1]
    cm[4*i+2,4*i+2]=xp[i+1]**2
    cm[4*i+2,4*i+3]=xp[i+1]**3
    cm[4*i+3,4*i]=0.0
    if i==len(xp)-2:
      cm[4*i+3,4*i+1]=0.0
      cm[4*i+3,4*i+2]=1.0
      cm[4*i+3,4*i+3]=3*xright
    else:
      cm[4*i+3,4*i+1]=1.0
      cm[4*i+3,4*i+2]=2*xp[i+1]
      cm[4*i+3,4*i+3]=3*xp[i+1]**2
      cm[4*i+3,4*i+4]=0.0
      cm[4*i+3,4*i+5]=-1.0
      cm[4*i+3,4*i+6]=-2*xp[i+1]
      cm[4*i+3,4*i+7]=-3*xp[i+1]**2
  cs=zeros(4*(len(xp)-1))
  for i in range(0,len(xp)-1):
    cs[4*i+1]=yp[i]
    cs[4*i+2]=yp[i+1]
  cc=solve(cm,cs)
  yq=zeros(len(xq))
  for q in range(0,len(xq)):
    i=len(xp)-2
    while (i>0) and (xq[q]<xp[i]):
      i-=1
    yq[q]=cc[4*i+3]*xq[q]**3+cc[4*i+2]*xq[q]**2+cc[4*i+1]*xq[q]+cc[4*i]
  return yq
% nicht im Command Window
function yq=mysplines(xp,yp,xleft,xright,xq)
cm=zeros(4*(length(xp)-1),4*(length(xp)-1));
for i=1:length(xp)-1;
cm(4*i-3,4*i-3)=0;
cm(4*i-3,4*i-2)=0;
cm(4*i-3,4*i-1)=1;
if i==1
cm(4*i-3,4*i)=3*xleft;
else
cm(4*i-3,4*i)=3*xp(i);
cm(4*i-3,4*i-7)=0;
cm(4*i-3,4*i-6)=0;
cm(4*i-3,4*i-5)=-1;
cm(4*i-3,4*i-4)=-3*xp(i);
end
cm(4*i-2,4*i-3)=1;
cm(4*i-2,4*i-2)=xp(i);
cm(4*i-2,4*i-1)=xp(i)^2;
cm(4*i-2,4*i)=xp(i)^3;
cm(4*i-1,4*i-3)=1;
cm(4*i-1,4*i-2)=xp(i+1);
cm(4*i-1,4*i-1)=xp(i+1)^2;
cm(4*i-1,4*i)=xp(i+1)^3;
cm(4*i,4*i-3)=0;
if i==length(xp)-1;
cm(4*i,4*i-2)=0;
cm(4*i,4*i-1)=1;
cm(4*i,4*i)=3*xright;
else
cm(4*i,4*i-2)=1;
cm(4*i,4*i-1)=2*xp(i+1);
cm(4*i,4*i)=3*xp(i+1)^2;
cm(4*i,4*i+1)=0;
cm(4*i,4*i+2)=-1;
cm(4*i,4*i+3)=-2*xp(i+1);
cm(4*i,4*i+4)=-3*xp(i+1)^2;
end
end
cs=zeros(4*(length(xp)-1));
for i=1:length(xp)-1
cs(4*i-2)=yp(i);
cs(4*i-1)=yp(i+1);
end
cc=linsolve(cm,cs);
yq=zeros(1,length(xq));
for q=1:length(xq)
i=length(xp)-1;
while ((i>1) && (xq(q)<xp(i)))
i=i-1;
end
yq(q)=cc(4*i)*xq(q)^3+cc(4*i-1)*xq(q)^2+cc(4*i-2)*xq(q)+cc(4*i-3);
end
Auswahl geeigneter Interpolationen in Abhängigkeit von der Anzahl der zu interpolierenden Extrempunkte
def interpol(mea,mev,mee):
  if len(mea)>2:
    #Splines, zweite Ableitung an den Raendern!! null
    meq=mysplines(mea,mev,0,N-1,arange(0,N))
  else:
    if len(mea)==2:
      #lineare Inter-/Extrapolation
      meq=arange(0,N)*(mev[1]-mev[0])/float(mea[1]-mea[0])+mev[0]-mea[0]*(mev[1]-mev[0])/float(mea[1]-mea[0])
    else:
      if len(mea)==1:
        #Konstante
        meq=ones(N)*mev[0]
      else:
        #Ersatzwert
        meq=ones(N)*mee
      
    
  return meq
% nicht im Command Window
function meq=interpol(N,mea,mev,mee)
if length(mea)>2
meq=mysplines(mea,mev,0,N-1,0:N-1);
else
if length(mea)==2
meq=(0:N-1)*(mev(2)-mev(1))/(mea(2)-mea(1))+mev(1)-mea(1)*(mev(2)-mev(1))/(mea(2)-mea(1));
else
if length(mea)==1
meq=ones(1,N)*mev(1);
else
meq=ones(1,N)*mee;
end
end
end
EMD, Zerlegung in intrinsische Moden
xm=copy(x)
modes=zeros([0,N])
for m in range(0,N):
  if max(xm)>min(xm):
    mmm=zeros(N)
    for l in range(0,100): #Anzahl der Siebschritte
      xmm=xm-mmm
      mxv=array([])
      mxa=array([])
      for i in range(1,N-1):
        if (xmm[i]>xmm[i-1]) and (xmm[i]>xmm[i+1]):
          mxv=append(mxv,xmm[i])
          mxa=append(mxa,i)
        
      mnv=array([])
      mna=array([])
      for i in range(1,N-1):
        if (xmm[i]<xmm[i-1]) and (xmm[i]<xmm[i+1]):
          mnv=append(mnv,xmm[i])
          mna=append(mna,i)
        
      mxq=interpol(mxa,mxv,max(xmm))
      mnq=interpol(mna,mnv,min(xmm))
      mmm+=0.5*(mxq+mnq)
    modes=append(modes,xmm.reshape((1,N)),axis=0)
    xm=mmm
  

plot(arange(0,N),modes[0,:],arange(0,N),modes[1,:],arange(0,N),modes[2,:])
show()
xm=x;
modes=zeros(0,N);
for m=1:N
if max(xm)>min(xm)
mmm=zeros(1,N);
for l=1:100 % Anzahl der Siebschritte
xmm=xm-mmm;
mxv=[];
mxa=[];
for i=2:N-1
if ((xmm(i)>xmm(i-1)) && (xmm(i)>xmm(i+1)))
mxv=[mxv,xmm(i)];
mxa=[mxa,i];
end
end
mnv=[];
mna=[];
for i=2:N-1
if ((xmm(i)<xmm(i-1)) && (xmm(i)<xmm(i+1)))
mnv=[mnv,xmm(i)];
mna=[mna,i];
end
end
% interpol(N,mxa,mxv,max(xmm))
if length(mxa)>2
cm=zeros(4*(length(mxa)-1),4*(length(mxa)-1));
for i=1:length(mxa)-1;
cm(4*i-3,4*i-3)=0;
cm(4*i-3,4*i-2)=0;
cm(4*i-3,4*i-1)=1;
if i==1
cm(4*i-3,4*i)=3*0;
else
cm(4*i-3,4*i)=3*mxa(i);
cm(4*i-3,4*i-7)=0;
cm(4*i-3,4*i-6)=0;
cm(4*i-3,4*i-5)=-1;
cm(4*i-3,4*i-4)=-3*mxa(i);
end
cm(4*i-2,4*i-3)=1;
cm(4*i-2,4*i-2)=mxa(i);
cm(4*i-2,4*i-1)=mxa(i)^2;
cm(4*i-2,4*i)=mxa(i)^3;
cm(4*i-1,4*i-3)=1;
cm(4*i-1,4*i-2)=mxa(i+1);
cm(4*i-1,4*i-1)=mxa(i+1)^2;
cm(4*i-1,4*i)=mxa(i+1)^3;
cm(4*i,4*i-3)=0;
if i==length(mxa)-1;
cm(4*i,4*i-2)=0;
cm(4*i,4*i-1)=1;
cm(4*i,4*i)=3*(N-1);
else
cm(4*i,4*i-2)=1;
cm(4*i,4*i-1)=2*mxa(i+1);
cm(4*i,4*i)=3*mxa(i+1)^2;
cm(4*i,4*i+1)=0;
cm(4*i,4*i+2)=-1;
cm(4*i,4*i+3)=-2*mxa(i+1);
cm(4*i,4*i+4)=-3*mxa(i+1)^2;
end
end
cs=zeros(4*(length(mxa)-1));
for i=1:length(mxa)-1
cs(4*i-2)=mxv(i);
cs(4*i-1)=mxv(i+1);
end
cc=linsolve(cm,cs);
yq=zeros(1,N);
for q=0:N-1
i=length(mxa)-1;
while ((i>1) && (q<mxa(i)))
i=i-1;
end
mxq(q+1)=cc(4*i)*q^3+cc(4*i-1)*q^2+cc(4*i-2)*q+cc(4*i-3);
end
else
if length(mxa)==2
mxq=(0:N-1)*(mxv(2)-mxv(1))/(mxa(2)-mxa(1))+mxv(1)-mxa(1)*(mxv(2)-mxv(1))/(mxa(2)-mxa(1));
else
if length(mxa)==1
mxq=ones(1,N)*mxv(1);
else
mxq=ones(1,N)*max(xmm);
end
end
end
% interpol(N,mna,mnv,min(xmm))
if length(mna)>2
cm=zeros(4*(length(mna)-1),4*(length(mna)-1));
for i=1:length(mna)-1;
cm(4*i-3,4*i-3)=0;
cm(4*i-3,4*i-2)=0;
cm(4*i-3,4*i-1)=1;
if i==1
cm(4*i-3,4*i)=3*0;
else
cm(4*i-3,4*i)=3*mna(i);
cm(4*i-3,4*i-7)=0;
cm(4*i-3,4*i-6)=0;
cm(4*i-3,4*i-5)=-1;
cm(4*i-3,4*i-4)=-3*mna(i);
end
cm(4*i-2,4*i-3)=1;
cm(4*i-2,4*i-2)=mna(i);
cm(4*i-2,4*i-1)=mna(i)^2;
cm(4*i-2,4*i)=mna(i)^3;
cm(4*i-1,4*i-3)=1;
cm(4*i-1,4*i-2)=mna(i+1);
cm(4*i-1,4*i-1)=mna(i+1)^2;
cm(4*i-1,4*i)=mna(i+1)^3;
cm(4*i,4*i-3)=0;
if i==length(mna)-1;
cm(4*i,4*i-2)=0;
cm(4*i,4*i-1)=1;
cm(4*i,4*i)=3*(N-1);
else
cm(4*i,4*i-2)=1;
cm(4*i,4*i-1)=2*mna(i+1);
cm(4*i,4*i)=3*mna(i+1)^2;
cm(4*i,4*i+1)=0;
cm(4*i,4*i+2)=-1;
cm(4*i,4*i+3)=-2*mna(i+1);
cm(4*i,4*i+4)=-3*mna(i+1)^2;
end
end
cs=zeros(4*(length(mna)-1));
for i=1:length(mna)-1
cs(4*i-2)=mnv(i);
cs(4*i-1)=mnv(i+1);
end
cc=linsolve(cm,cs);
yq=zeros(1,N);
for q=0:N-1
i=length(mna)-1;
while ((i>1) && (q<mna(i)))
i=i-1;
end
mnq(q+1)=cc(4*i)*q^3+cc(4*i-1)*q^2+cc(4*i-2)*q+cc(4*i-3);
end
else
if length(mna)==2
mnq=(0:N-1)*(mnv(2)-mnv(1))/(mna(2)-mna(1))+mnv(1)-mna(1)*(mnv(2)-mnv(1))/(mna(2)-mna(1));
else
if length(mna)==1
mnq=ones(1,N)*mnv(1);
else
mnq=ones(1,N)*min(xmm);
end
end
end
mmm=mmm+0.5*(mxq+mnq);
end
modes=[modes;xmm];
xm=mmm;
end
end
plot(0:N-1,modes(1,:),0:N-1,modes(2,:),0:N-1,modes(3,:))

Proper Orthogonal Decomposition (POD)

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
Generieren eines Beispieldatensatzes von 100 Snapshots mit 40 Merkmalen (Summe einer harmonischen Schwingung und einer harmonischen Welle mit Rauschen)
N=100
p=40
x=zeros([N,p])
sigma1=15.0
sigma2=15.0
for i in range(0,N):
  for j in range(0,p):
    t1=(j+0.5-p/2.0)+0.3*(i+0.5-N/2.0)
    t2=(j+0.5-p/2.0)+0.0*(i+0.5-N/2.0)
    x[i,j]=+0.5*sin(0.6*t1)*exp(-(j+0.5-p/2.0)**2/(2*sigma1**2))+1.0*sin(0.45*t2)*exp(-(j+0.5-p/2.0)**2/(2*sigma2**2))+0.1*standard_normal()
  

imshow(x,cmap='gray')
show()
N=100;
p=40;
x=zeros(N,p);
sigma1=15;
sigma2=15;
for i=1:N
for j=1:p
t1=(j-0.5-p/2)+0.3*(i-0.5-N/2);
t2=(j-0.5-p/2)+0.0*(i-0.5-N/2);
x(i,j)=+0.5*sin(0.6*t1)*exp(-(j-0.5-p/2)^2/(2*sigma1^2))+1.0*sin(0.45*t2)*exp(-(j-0.5-p/2)^2/(2*sigma2^2))+0.1*randn();
end
end
imshow(x,[min(min(x)),max(max(x))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Korrelationsmatrix
C=dot(x.T,x)
imshow(C+amin(C),cmap='gray')
show()
C=x'*x;
imshow(C,[min(min(C)),max(max(C))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Eigenwerte und Eigenvektoren (Moden)
from numpy.linalg import *
eval,evek=eig(C)
evek=real(evek)
for i in range(0,p-1):
  for j in range(i+1,p):
    if abs(eval[i])<abs(eval[j]):
      evalx=copy(eval[i])
      eval[i]=eval[j]
      eval[j]=evalx
      evekx=copy(evek[:,i])
      evek[:,i]=evek[:,j]
      evek[:,j]=evekx
    
  

plot(arange(0,p),eval,arange(0,p),eval,'o')
show()
imshow(evek+amin(evek),cmap='gray')
show()
[evek,eval]=eig(C);
evek=real(evek);
for i=1:p-1
for j=i+1:p
if (abs(eval(i,i))<abs(eval(j,j)))
evalx=eval(i,i);
eval(i,i)=eval(j,j);
eval(j,j)=evalx;
evekx=evek(:,i);
evek(:,i)=evek(:,j);
evek(:,j)=evekx;
end
end
end
plot(1:p,diag(eval),1:p,diag(eval),'o')
pause
imshow(evek,[min(min(evek)),max(max(evek))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Transformation
coef=zeros([N,p])
for i in range(0,N):
  for j in range(0,p):
    coef[i,j]=dot(x[i,:],evek[:,j]);
  

imshow(coef,cmap='gray')
show()
coef=zeros(N,p);
for i=1:N
for j=1:p
coef(i,j)=x(i,:)*evek(:,j);
end
end
imshow(coef,[min(min(coef)),max(max(coef))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Rekonstruktion eines bestimmten Modes
xrek=zeros([N,p])
m=0 # Index null für ersten Mode
for i in range(0,N):
  for j in range(0,p):
    xrek[i,j]=dot(coef[i,m],evek[j,m])
  

imshow(xrek,cmap='gray')
show()
xrek=zeros(N,p);
m=1; % Index eins für ersten Mode
for i=1:N
for j=1:p
xrek(i,j)=coef(i,m)*evek(j,m)';
end
end
imshow(xrek,[min(min(xrek)),max(max(xrek))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
vollständige Rekonstruktion
xrek=zeros([N,p])
for i in range(0,N):
  for j in range(0,p):
    xrek[i,j]=dot(coef[i,:],evek[j,:])
  

imshow(xrek,cmap='gray')
show()
xrek=zeros(N,p);
for i=1:N
for j=1:p
xrek(i,j)=coef(i,:)*evek(j,:)';
end
end
imshow(xrek,[min(min(xrek)),max(max(xrek))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Rekonstruktion aus den ersten drei Moden
xrek=zeros([N,p])
m=range(0,3) # für die ersten drei Moden
for i in range(0,N):
  for j in range(0,p):
    xrek[i,j]=dot(coef[i,m],evek[j,m])
  

imshow(xrek,cmap='gray')
show()
xrek=zeros(N,p);
m=[1:3]; % für die ersten drei Moden
for i=1:N
for j=1:p
xrek(i,j)=coef(i,m)*evek(j,m)';
end
end
imshow(xrek,[min(min(xrek)),max(max(xrek))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])

Snapshot-POD

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
Generieren eines Beispieldatensatzes von 100 Snapshots mit 40 Merkmalen (Summe einer harmonischen Schwingung und einer harmonischen Welle mit Rauschen)
N=100
p=40
x=zeros([p,N])
sigma1=1000.0
sigma2=1000.0
for i in range(0,p):
  for j in range(0,N):
    t1=1.0*(i+0.5-p/2.0)+0.3*(j+0.5-N/2.0)
    t2=0.0*(i+0.5-p/2.0)-0.5*(j+0.5-N/2.0)
    x[i,j]=+0.2*sin(0.6*t1)*exp(-(j+0.5-N/2.0)**2/(2*sigma1**2))+0.6*sin(0.45*t2)*cos(pi*(i+0.5-p/2.0)/float(p))*exp(-(j+0.5-N/2.0)**2/(2*sigma2**2))+0.01*standard_normal()
  

imshow(x,cmap='gray')
show()
N=100;
p=40;
x=zeros(p,N);
sigma1=1000;
sigma2=1000;
for i=1:p
for j=1:N
t1=1.0*(i-0.5-p/2)+0.3*(j-0.5-N/2);
t2=0.0*(i-0.5-p/2)-0.5*(j-0.5-N/2);
x(i,j)=+0.2*sin(0.6*t1)*exp(-(j-0.5-N/2)^2/(2*sigma1^2))+0.6*sin(0.45*t2)*cos(pi*(i-0.5-p/2)/p)*exp(-(j-0.5-N/2)^2/(2*sigma2^2))+0.01*randn();
end
end
imshow(x,[min(min(x)),max(max(x))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Korrelationsmatrix
C=dot(x.T,x)
imshow(C+amin(C),cmap='gray')
show()
C=x'*x;
imshow(C,[min(min(C)),max(max(C))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Eigenwerte und Eigenvektoren (Moden)
from numpy.linalg import *
eval,evek=eig(C)
evek=real(evek)
for i in range(0,N-1):
  for j in range(i+1,N):
    if abs(eval[i])<abs(eval[j]):
      evalx=copy(eval[i])
      eval[i]=eval[j]
      eval[j]=evalx
      evekx=copy(evek[:,i])
      evek[:,i]=evek[:,j]
      evek[:,j]=evekx
    
  

plot(arange(0,N),eval,arange(0,N),eval,'o')
show()
imshow(evek+amin(evek),cmap='gray')
show()
[evek,eval]=eig(C);
evek=real(evek);
for i=1:N-1
for j=i+1:N
if (abs(eval(i,i))<abs(eval(j,j)))
evalx=eval(i,i);
eval(i,i)=eval(j,j);
eval(j,j)=evalx;
evekx=evek(:,i);
evek(:,i)=evek(:,j);
evek(:,j)=evekx;
end
end
end
plot(1:N,diag(eval),1:N,diag(eval),'o')
pause
imshow(evek,[min(min(evek)),max(max(evek))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Transformation
coef=zeros([p,N])
for i in range(0,p):
  for j in range(0,N):
    coef[i,j]=dot(x[i,:],evek[:,j]);
  

imshow(coef,cmap='gray')
show()
coef=zeros(p,N);
for i=1:p
for j=1:N
coef(i,j)=x(i,:)*evek(:,j);
end
end
imshow(coef,[min(min(coef)),max(max(coef))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Rekonstruktion eines bestimmten Modes
xrek=zeros([p,N])
m=0 # Index null für ersten Mode
for i in range(0,p):
  for j in range(0,N):
    xrek[i,j]=dot(coef[i,m],evek[j,m])
  

imshow(xrek,cmap='gray')
show()
xrek=zeros(p,N);
m=1; % Index eins für ersten Mode
for i=1:p
for j=1:N
xrek(i,j)=coef(i,m)*evek(j,m)';
end
end
imshow(xrek,[min(min(xrek)),max(max(xrek))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
vollständige Rekonstruktion
xrek=zeros([p,N])
for i in range(0,p):
  for j in range(0,N):
    xrek[i,j]=dot(coef[i,:],evek[j,:])
  

imshow(xrek,cmap='gray')
show()
xrek=zeros(p,N);
for i=1:p
for j=1:N
xrek(i,j)=coef(i,:)*evek(j,:)';
end
end
imshow(xrek,[min(min(xrek)),max(max(xrek))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Rekonstruktion aus den ersten drei Moden
xrek=zeros([p,N])
m=range(0,3) # für die ersten drei Moden
for i in range(0,p):
  for j in range(0,N):
    xrek[i,j]=dot(coef[i,m],evek[j,m])
  

imshow(xrek,cmap='gray')
show()
xrek=zeros(p,N);
m=[1:3]; % für die ersten drei Moden
for i=1:p
for j=1:N
xrek(i,j)=coef(i,m)*evek(j,m)';
end
end
imshow(xrek,[min(min(xrek)),max(max(xrek))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])

komplex erweiterte Snapshot-POD

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
Generieren eines Beispieldatensatzes von 100 Snapshots mit 40 Merkmalen (Summe einer harmonischen Schwingung und einer harmonischen Welle mit Rauschen)
N=100
p=40
x=zeros([p,N])
sigma1=1000.0
sigma2=1000.0
for i in range(0,p):
  for j in range(0,N):
    t1=1.0*(i+0.5-p/2.0)+0.3*(j+0.5-N/2.0)
    t2=0.0*(i+0.5-p/2.0)-0.5*(j+0.5-N/2.0)
    x[i,j]=+0.2*sin(0.6*t1)*exp(-(j+0.5-N/2.0)**2/(2*sigma1**2))+0.6*sin(0.45*t2)*cos(pi*(i+0.5-p/2.0)/float(p))*exp(-(j+0.5-N/2.0)**2/(2*sigma2**2))+0.01*standard_normal()
  

imshow(x,cmap='gray')
show()
N=100;
p=40;
x=zeros(p,N);
sigma1=1000;
sigma2=1000;
for i=1:p
for j=1:N
t1=1.0*(i-0.5-p/2)+0.3*(j-0.5-N/2);
t2=0.0*(i-0.5-p/2)-0.5*(j-0.5-N/2);
x(i,j)=+0.2*sin(0.6*t1)*exp(-(j-0.5-N/2)^2/(2*sigma1^2))+0.6*sin(0.45*t2)*cos(pi*(i-0.5-p/2)/p)*exp(-(j-0.5-N/2)^2/(2*sigma2^2))+0.01*randn();
end
end
imshow(x,[min(min(x)),max(max(x))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Komplexe Erweiterung
from scipy.signal import *
xc=zeros([p,N])+0j
for i in range(0,p):
  xc[i]=hilbert(x[i]-mean(x[i]))+mean(x[i])

imshow(imag(xc),cmap='gray')
show()
xc=zeros(p,N)+0j;
for i=1:p
Xi=fft(x(i,:));
Xi(1)=0;
for j=2:fix((N+1)/2)
Xi(j)=-1j*Xi(j);
end
for j=fix((N+3)/2):fix((N+2)/2)
Xi(j)=0;
end
for j=fix((N+4)/2):N
Xi(j)=+1j*Xi(j);
end
xc(i,:)=x(i,:)+1j*ifft(Xi);
end
imshow(imag(xc),[min(min(imag(xc))),max(max(imag(xc)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Korrelationsmatrix
C=dot(xc.T,conj(xc))
imshow(real(C)+amin(real(C)),cmap='gray')
show()
imshow(imag(C)+amin(imag(C)),cmap='gray')
show()
C=conj(xc')*conj(xc); % Achtung! Matlab/Octave erzeugt bei komplexem xc durch xc' die konjuguert komplexe Transponierte.
imshow(real(C),[min(min(real(C))),max(max(real(C)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
pause
imshow(imag(C),[min(min(imag(C))),max(max(imag(C)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Eigenwerte und Eigenvektoren (Moden)
from numpy.linalg import *
eval,evek=eig(C)
for i in range(0,N-1):
  for j in range(i+1,N):
    if abs(eval[i])<abs(eval[j]):
      evalx=copy(eval[i])
      eval[i]=eval[j]
      eval[j]=evalx
      evekx=copy(evek[:,i])
      evek[:,i]=evek[:,j]
      evek[:,j]=evekx
    
  

plot(arange(0,N),eval,arange(0,N),eval,'o')
show()
imshow(real(evek)+amin(real(evek)),cmap='gray')
show()
imshow(imag(evek)+amin(imag(evek)),cmap='gray')
show()
[evek,eval]=eig(C);
for i=1:N-1
for j=i+1:N
if (abs(eval(i,i))<abs(eval(j,j)))
evalx=eval(i,i);
eval(i,i)=eval(j,j);
eval(j,j)=evalx;
evekx=evek(:,i);
evek(:,i)=evek(:,j);
evek(:,j)=evekx;
end
end
end
plot(1:N,diag(eval),1:N,diag(eval),'o')
pause
imshow(real(evek),[min(min(real(evek))),max(max(real(evek)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
pause
imshow(imag(evek),[min(min(imag(evek))),max(max(imag(evek)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Transformation
coef=zeros([p,N])+0j
for i in range(0,p):
  for j in range(0,N):
    coef[i,j]=dot(xc[i,:],conj(evek[:,j]));
  

imshow(real(coef),cmap='gray')
show()
imshow(imag(coef),cmap='gray')
show()
coef=zeros(p,N)+0j;
for i=1:p
for j=1:N
coef(i,j)=xc(i,:)*conj(evek(:,j));
end
end
imshow(real(coef),[min(min(real(coef))),max(max(real(coef)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
pause
imshow(imag(coef),[min(min(imag(coef))),max(max(imag(coef)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Rekonstruktion eines bestimmten Modes
xrek=zeros([p,N])+0j
m=0 # Index null für ersten Mode
for i in range(0,p):
  for j in range(0,N):
    xrek[i,j]=dot(coef[i,m],evek[j,m])
  

imshow(real(xrek),cmap='gray')
show()
imshow(imag(xrek),cmap='gray')
show()
xrek=zeros(p,N)+0j;
m=1; % Index eins für ersten Mode
for i=1:p
for j=1:N
xrek(i,j)=coef(i,m)*conj(evek(j,m)'); % Achtung! Matlab/Octave erzeugt bei komplexem xc durch xc' die konjuguert komplexe Transponierte.
end
end
imshow(real(xrek),[min(min(real(xrek))),max(max(real(xrek)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
pause
imshow(imag(xrek),[min(min(imag(xrek))),max(max(imag(xrek)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
vollständige Rekonstruktion
xrek=zeros([p,N])+0j
for i in range(0,p):
  for j in range(0,N):
    xrek[i,j]=dot(coef[i,:],evek[j,:])
  

imshow(real(xrek),cmap='gray')
show()
imshow(imag(xrek),cmap='gray')
show()
xrek=zeros(p,N)+0j;
for i=1:p
for j=1:N
xrek(i,j)=coef(i,:)*conj(evek(j,:)'); % Achtung! Matlab/Octave erzeugt bei komplexem xc durch xc' die konjuguert komplexe Transponierte.
end
end
imshow(real(xrek),[min(min(real(xrek))),max(max(real(xrek)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
pause
imshow(imag(xrek),[min(min(imag(xrek))),max(max(imag(xrek)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Rekonstruktion aus den ersten drei Moden
xrek=zeros([p,N])+0j
m=range(0,3) # für die ersten drei Moden
for i in range(0,p):
  for j in range(0,N):
    xrek[i,j]=dot(coef[i,m],evek[j,m])
  

imshow(real(xrek),cmap='gray')
show()
imshow(imag(xrek),cmap='gray')
show()
xrek=zeros(p,N)+0j;
m=[1:3]; % für die ersten drei Moden
for i=1:p
for j=1:N
xrek(i,j)=coef(i,m)*conj(evek(j,m)'); % Achtung! Matlab/Octave erzeugt bei komplexem xc durch xc' die konjuguert komplexe Transponierte.
end
end
imshow(real(xrek),[min(min(real(xrek))),max(max(real(xrek)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
pause
imshow(imag(xrek),[min(min(imag(xrek))),max(max(imag(xrek)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])