Signal- und Messdatenverarbeitung


Testcodes für Kapitel 16: Wavelets

systematische Zusammenstellung von kombinierbaren Programmabschnitten zur Signal- und Messdatenverarbeitung

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
#python -m pip install --upgrade PyWavelets
#python -m pip install --upgrade numpy
#python -m pip install --upgrade matplotlib
#python -m pip install --upgrade scipy
#python -m pip install --upgrade PIL
#pip install --upgrade PyWavelets
#pip install --upgrade numpy
#pip install --upgrade matplotlib
#pip install --upgrade scipy
#pip install --upgrade PIL
from numpy import *
from matplotlib.pyplot import *
%Matlab: erfordert die Wavelet Toolbox
Signal
from sys import *
if version_info.major==2:
  from urllib2 import *

if version_info.major==3:
  from urllib import *

tx=array(loadtxt(urlopen('http://nambis.bplaced.de/nambis/SMDV/echolot.dat'),delimiter='\t'))
t=tx[:,0]
x=tx[:,1]
plot(t,x)
show()
tx=textscan(urlread('http://nambis.bplaced.de/nambis/SMDV/echolot.dat'),'%f\t%f','EndOfLine','\n');
t=tx{1};
x=tx{2};
plot(t,x)
kontinuierliche WT
from scipy.signal import *
W=array(cwt(x,ricker,arange(1,31))) #ricker=mexican hat
#W=cwt(x,morlet,arange(1,31))
#W=cwt(x,daub,arange(1,31))
#W=cwt(x,qmf,arange(1,31))
#W=cwt(x,cascade,arange(1,31))
#W=cwt(x,cwt,arange(1,31))
imshow(W,extent=[-1,1,31,1],cmap='PRGn',aspect='auto',vmax=abs(W).max(),vmin=-abs(W).max())
show()
for i in range(0,30):
  plot(W[i])
  show()

Matlab:
W=cwt(x,'morse'); % Voreinstellung
%W=cwt(x,'amor'); % Morlet
%W=cwt(x,'bump');
imshow(real(W),[min(min(real(W))),max(max(real(W)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
waitforbuttonpress
imshow(imag(W),[min(min(imag(W))),max(max(imag(W)))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
waitforbuttonpress
for i=1:2:length(W(:,1))
plot([0:length(W(i,:))-1],real(W(i,:)),[0:length(W(i,:))-1],imag(W(i,:)))
waitforbuttonpress
end
Octave:
pkg load ltfat
W=fwt(x,'db8',10);
weiter???
diskrete WT (eindimensional, Datenreduktion, Rekonstruierbarkeit)
from pywt import *
cA,cD=dwt(x,'db2')
xr=idwt(cA,cD,'db2')
xa=idwt(cA,zeros(size(cD)),'db2')
plot(t,x,t,xr,t,xa)
show()
cA,cD=dwt(x,'haar')
xr=idwt(cA,cD,'haar')
xa=idwt(cA,zeros(size(cD)),'haar')
plot(t,x,t,xr,t,xa)
show()
Matlab:
[cA,cD]=dwt(x,'db2');
xr=idwt(cA,cD,'db2');
xa=idwt(cA,zeros(size(cD)),'db2');
plot(t,x,t,xr,t,xa)
waitforbuttonpress
[cA,cD]=dwt(x,'haar');
xr=idwt(cA,cD,'haar');
xa=idwt(cA,zeros(size(cD)),'haar');
plot(t,x,t,xr,t,xa)
Octave:
fehlt noch
Bild
from PIL.Image import *
b=open(urlopen('http://nambis.bplaced.de/nambis/SMDV/schneemann.png'))
type(b)
imshow(b,cmap='gray')
show()
b=imread('http://nambis.bplaced.de/nambis/SMDV/schneemann.png');
imshow(b)
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
diskrete WT (zweidimensional, Datenreduktion, Rekonstruierbarkeit)
from pywt import *
coeffs=wavedec2(b,'haar')
imshow(coeffs[0],cmap='gray')
show()
for i in range(1,len(coeffs)):
  imshow(coeffs[i][2],cmap='gray')
  show()

br=waverec2(coeffs,'haar')
imshow(br,cmap='gray')
show()
for i in range(6,len(coeffs)): # Datenreduktion, Ebenen 6 und 7 entfernt
  for j in range(0,len(coeffs[i])):
    for k in range(0,len(coeffs[i][j])):
     for l in range(0,len(coeffs[i][j][k])):
        coeffs[i][j][k,l]=0
      
    
  

br=waverec2(coeffs,'haar')
imshow(br,cmap='gray')
show()
Matlab:
[coeffs,info]=wavedec2(b,7,'haar');
br=waverec2(coeffs,info,'haar');
imshow(br,[min(min(br)),max(max(br))])
imgps=get(gcf,'Position');
set(gcf,'Position',[imgps(1),imgps(2),560,420])
Octave:
fehlt noch

dyadische Transformation

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
#python -m pip install --upgrade numpy
#python -m pip install --upgrade matplotlib
#python -m pip install --upgrade PIL
#python -m pip install --upgrade scipy
#pip install --upgrade numpy
#pip install --upgrade matplotlib
#pip install --upgrade PIL
#pip install --upgrade scipy
from numpy import *
from matplotlib.pyplot import *
from PIL.Image import *
Originalbild laden und anzeigen
from sys import *
if version_info.major==2:
  from urllib2 import *

if version_info.major==3:
  from urllib import *

borig=open(urlopen('http://nambis.bplaced.de/nambis/SMDV/schneemann.png'))
imshow(borig,cmap='gray')
show()
fehlt noch
dyadische Zerlegung mit Wavelets
from pywt import *
[LL1,[HL1,LH1,HH1]]=dwt2(borig,'haar')
bcoef=concatenate((concatenate((LL1,LH1),axis=1),concatenate((HL1,HH1),axis=1)),axis=0)
imshow(bcoef, cmap='gray')
show()
fehlt noch
vollständige Rekonstruktion des Bildes
brec=idwt2([LL1,[HL1,LH1,HH1]],'haar')
imshow(brec,cmap='gray')
show()
fehlt noch
dyadische Zerlegung mit der DFT
from numpy.fft import *
w,h=borig.size
Borig=fft2(borig)
LL1=concatenate((concatenate((Borig[0:h//4,0:w//4],Borig[h-(h//4):h,0:w//4]),axis=0),concatenate((Borig[0:h//4,w-(w//4):w],Borig[h-(h//4):h,w-(w//4):w]),axis=0)),axis=1)
LH1=concatenate((Borig[0:h//4,w//4:w-(w//4)],Borig[h-(h//4):h,w//4:w-(w//4)]),axis=0)
HL1=concatenate((Borig[h//4:h-(h//4),0:w//4],Borig[h//4:h-(h//4),w-(w//4):w]),axis=1)
HH1=Borig[h//4:h-(h//4),w//4:w-(w//4)]
ll1=real(ifft2(LL1))
lh1=real(ifft2(LH1))
hl1=real(ifft2(HL1))
hh1=real(ifft2(HH1))
bcoef=concatenate((concatenate((ll1,lh1),axis=1),concatenate((hl1,hh1),axis=1)),axis=0)
imshow(bcoef,cmap='gray')
show()
fehlt noch
vollständige Rekonstruktion des Bildes
LL1=fft2(bcoef[0:h//2,0:w//2])
LH1=fft2(bcoef[0:h//2,w//2:w])
HL1=fft2(bcoef[h//2:h,0:w//2])
HH1=fft2(bcoef[h//2:h,w//2:w])
Brec=concatenate((concatenate((LL1[0:h//4,0:w//4],LH1[0:h//4,:],LL1[0:h//4,w//2-(w//4):w//2]),axis=1),concatenate((HL1[:,0:w//4],HH1[:,:],HL1[:,w//2-(w//4):w//2]),axis=1),concatenate((LL1[h//2-(h//4):h//2,0:w//4],LH1[h//2-(h//4):h//2,:],LL1[h//2-(h//4):h//2,w//2-(w//4):w//2]),axis=1)),axis=0)
brec=real(ifft2(Brec))
imshow(brec,cmap='gray')
show()
fehlt noch
dyadische Zerlegung mit der DCT
from scipy.fftpack import *
w,h=borig.size
Borig=dct(dct(borig,axis=1),axis=0)
LL1=Borig[0:h//2,0:w//2]
LH1=Borig[0:h//2,w-(w//2):w]
HL1=Borig[h-(h//2):h,0:w//2]
HH1=Borig[h-(h//2):h,w-(w//2):w]
ll1=idct(idct(LL1,axis=0),axis=1)
lh1=idct(idct(LH1,axis=0),axis=1)
hl1=idct(idct(HL1,axis=0),axis=1)
hh1=idct(idct(HH1,axis=0),axis=1)
bcoef=concatenate((concatenate((ll1,lh1),axis=1),concatenate((hl1,hh1),axis=1)),axis=0)
imshow(bcoef,cmap='gray')
show()
fehlt noch
vollständige Rekonstruktion des Bildes
LL1=dct(dct(bcoef[0:h//2,0:w//2],axis=1),axis=0)
LH1=dct(dct(bcoef[0:h//2,w//2:w],axis=1),axis=0)
HL1=dct(dct(bcoef[h//2:h,0:w//2],axis=1),axis=0)
HH1=dct(dct(bcoef[h//2:h,w//2:w],axis=1),axis=0)
Brec=concatenate((concatenate((LL1,LH1),axis=1),concatenate((HL1,HH1),axis=1)),axis=0)
brec=idct(idct(Brec,axis=0),axis=1)
imshow(brec,cmap='gray')
show()
fehlt noch