Signal- und Messdatenverarbeitung

Testcodes für Kapitel 20: Quellenanalyse

systematische Zusammenstellung von kombinierbaren Programmabschnitten zur Signal- und Messdatenverarbeitung

Empirical Mode Decomposition (EMD)

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
#python -m pip install --upgrade numpy
#python -m pip install --upgrade matplotlib
#pip install --upgrade numpy
#pip install --upgrade matplotlib
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
from numpy.linalg import *
Generieren eines Signals
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(coef,yp,xleft,xrekight,xq):
  cm=zeros([4*(len(coef)-1),4*(len(coef)-1)])
  for i in range(0,len(coef)-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*coef[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*coef[i]
    cm[4*i+1,4*i]=1.0
    cm[4*i+1,4*i+1]=coef[i]
    cm[4*i+1,4*i+2]=coef[i]**2
    cm[4*i+1,4*i+3]=coef[i]**3
    cm[4*i+2,4*i]=1.0
    cm[4*i+2,4*i+1]=coef[i+1]
    cm[4*i+2,4*i+2]=coef[i+1]**2
    cm[4*i+2,4*i+3]=coef[i+1]**3
    cm[4*i+3,4*i]=0.0
    if i==len(coef)-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*xrekight
    else:
      cm[4*i+3,4*i+1]=1.0
      cm[4*i+3,4*i+2]=2*coef[i+1]
      cm[4*i+3,4*i+3]=3*coef[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*coef[i+1]
      cm[4*i+3,4*i+7]=-3*coef[i+1]**2
  cs=zeros(4*(len(coef)-1))
  for i in range(0,len(coef)-1):
    cs[4*i+1]=yp[i]
    cs[4*i+2]=yp[i+1]
  cc=solve(cm,cs)
  yq=zeros(size(xq))
  for q in range(0,len(xq)):
    i=len(coef)-2
    while (i>0) and (xq[q]<coef[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
%Matlab: 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
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
%Matlab: 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
end
EMD und Anzeige der gefundenen intrinsischen Moden
xm=copy(x)
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)
    plot(arange(0,N),xmm)
    show()
    xm=mmm
  
xm=x;
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
plot([0:N-1],xmm)
pause
%
xm=mmm;
end
end

Proper Orthogonal Decomposition (POD)

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
#python -m pip install --upgrade numpy
#python -m pip install --upgrade matplotlib
#pip install --upgrade numpy
#pip install --upgrade matplotlib
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
Generieren eines Feldes von 100 Vektoren mit 40 Merkmalen (Summe harmonischer Schwingungen)
M=100
p=40
x=zeros([M,p])
sigma1=15.0
sigma2=15.0
for i in range(0,M):
  for j in range(0,p):
    t1=(j+0.5-p/2.0)+0.3*(i+0.5-M/2.0)
    t2=(j+0.5-p/2.0)-0.0*(i+0.5-M/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()
M=100;
p=40;
x=zeros(M,p);
sigma1=15;
sigma2=15;
for i=1:M
for j=1:p
t1=(j-0.5-p/2)+0.3*(i-0.5-M/2);
t2=(j-0.5-p/2)+0.0*(i-0.5-M/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])
Kovarianzmatrix
C=zeros([p,p])
for i in range(0,p):
  for j in range(0,p):
    C[i,j]=sum(x[:,i]*x[:,j])
  

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([M,p])
for i in range(0,M):
  for j in range(0,p):
    coef[i,j]=dot(x[i,:],evek[:,j]);
  

imshow(coef,cmap='gray')
show()
coef=zeros(M,p);
for i=1:M
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])
vollständige Rekonstruktion
xrek=zeros([M,p])
for i in range(0,M):
  for j in range(0,p):
    xrek[i,j]=dot(coef[i,:],evek[j,:])
  

imshow(xrek,cmap='gray')
show()
xrek=zeros(M,p);
for i=1:M
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([M,p])
m=range(0,3) # für die ersten drei Moden
for i in range(0,M):
  for j in range(0,p):
    xrek[i,j]=dot(coef[i,m],evek[j,m])
  

imshow(xrek,cmap='gray')
show()
xrek=zeros(M,p);
m=[1:3]; % für die ersten drei Moden
for i=1:M
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])