Skip to content
Snippets Groups Projects
Commit 26ba52ae authored by Adam Friedrich Schrey's avatar Adam Friedrich Schrey
Browse files

Replace Matlab-Code with Python-Code and other smaller changes in V02_die_z-transformation.ipynb

parent d2858401
Branches
No related tags found
1 merge request!1merge develop into master
Show whitespace changes
Inline Side-by-side
notebooks/V02_die_z-transformation.ipynb
+ 595
− 200
View file @ 26ba52ae
%% Cell type:markdown id: tags:
# <span style='color:OrangeRed'>V2 - Die z-Transformation</span>
%% Cell type:code id: tags:
``` python
#%matplotlib widget
%matplotlib inline
import control
import ipywidgets as widgets
from IPython.display import display
from scipy import signal
import itertools
import matplotlib.pyplot as plt
import numpy as np
pi = np.pi
sin = np.sin
cos = np.cos
sqrt = np.sqrt
exp = np.exp
atan2 = np.arctan2
log10 = np.log10
def deconvolve(N,D,i):
f,R = signal.deconvolve(N,D)
print("deconvolve (im ersten Schritt N/D sonst R/D):")
print("f"+i+":"+str(f))
print("R:"+str(R))
return(f,R)
def convolve(R):
R = signal.convolve(R,[1, 0])
print("convolve with z^(-1): R=R*z^(-1):"+str(R)+"\n")
return(R)
def cont2dis_to_list(con,Ts,method):
dis = list(signal.cont2discrete(con,Ts,method))
dis[0] = dis[0].tolist()
dis[1] = dis[1].tolist()
for i in range(len(dis[0])):
if type(dis[0][i]) == list or type(dis[0][i]) == np.ndarray:
dis[0][i] = sum(dis[0][i])
for i in range(len(dis[1])):
if type(dis[1][i]) == list or type(dis[1][i]) == np.ndarray:
dis[1][i] = sum(dis[0][i])
return dis
def print_list_as_polynom(l,var):
n=0
for i in l:
if type(i) == list or type(i) == np.ndarray:
i = sum(i)
if n == 0:
print(" ",end="")
else:
print("+",end="")
print("("+str(np.around(i,3))+"*"+var+"^"+str(len(l)-n-1)+")", end = '')
n=n+1
print("")
def print_lists_as_division(l1,l2,var):
print_list_as_polynom(l1,var)
if len(l1) > len(l2):
n=len(l1)
else:
n=len(l2)
for i in range(10*n):
print("—",end="")
print("")
print_list_as_polynom(l2,var)
```
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #1 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Gegeben ist die Funktion im Zähler mit: $8z^{-1}$
Des Weiteren ist die Nennerfunktion gegeben mit:$1-3z^{-1}+2z^{-2}$
%% Cell type:code id: tags:
``` octave
N = [8 0]
D = [1 -3 2]
``` python
N = [8, 0]
D = [1, -3, 2]
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Da <code>N</code> und <code>D</code> Vektorpolynome sind , entspricht eine Entfaltung(.eng.: deconvolution) einer einfachen
Division. Die Matlab Funktion hierfür lautet:
%% Cell type:code id: tags:
``` octave
[fo R] = deconv(N,D)
``` python
f1,R = deconvolve(N,D,"1")
```
%% Output
deconvolve (im ersten Schritt N/D sonst R/D):
f1:[]
R:[8 0]
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Bei erstmaliger Anwendung der Funktion erhalten wir den ersten Koeffizienten der Folge <code>fo</code> und den
Rest ( eng.: remainder) <code>R</code>.
Der Rest <code>R</code> wird mit $z^{-1}$ multipliziert und weiter dividiert
%% Cell type:code id: tags:
``` octave
R = conv(R,[1 0])
[f2 R] = deconv(R,D)
``` python
R = convolve(R)
f2, R = deconvolve(R,D,"2")
```
%% Output
convolve with z^(-1): R=R*z^(-1):[8 0 0]
deconvolve (im ersten Schritt N/D sonst R/D):
f2:[8.]
R:[ 0. 24. -16.]
%% Cell type:code id: tags:
``` octave
R = conv(R,[1 0])
[f3 R] = deconv(R,D)
``` python
R = convolve(R)
f3, R = deconvolve(R,D,"3")
```
%% Output
convolve with z^(-1): R=R*z^(-1):[ 0. 24. -16. 0.]
deconvolve (im ersten Schritt N/D sonst R/D):
f3:[ 0. 24.]
R:[ 0. 0. 56. -48.]
%% Cell type:markdown id: tags:
### Die Potenzreihenentwicklung
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
<br><br>
<b> Mathematische Berechnungen: </b>
<br><br>
Es ist $f(k)$ für $k=0, 1, 2,...$ zu bestimmen, wenn
$F_z(z)= \frac{8z^{-1}}{(z-1)(z-2)}$
egeben ist.
<span style='color:OrangeRed'> Lösung:</span>
Durch Umschreiben folgt:
$F_z(z)= \frac{8z^{-1}}{1-3z^{-1}+2z^{-2}}$
Durch Division erhält man:
$F_z(z)=8z^{-1}+24z^{-2}+56z^{-3}+120z^{-4}+...$
Rechenweg (schriftliche Division):
![title](bilder/v02_rechenweg_division.png)
Aus dieser Beziehung können direkt die Werte der Zahlenfolge $f(k)$ abgelesen werden:
$f(0)=0$, $f(1)=8$, $f(2)=24$, $f(3)=56$, $f(4)=120$,...
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #2 </span>
%% Cell type:code id: tags:
%% Cell type:markdown id: tags:
``` octave
N = [8 2 3 2 1];
D = [1 0 0 0 0];
```
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
$F_z(z)= \frac{8z^{-4}+2z^{-3}+3z^{-2}+2z^{-1}+1}{z^{-4}}$
<h1>???</h1>
</div>
%% Cell type:code id: tags:
``` octave
[fo R] = deconv(N,D)
``` python
N = [8, 2, 3, 2, 1]
D = [1, 0, 0, 0, 0]
```
%% Cell type:code id: tags:
``` octave
R = conv(R,[1 0])
[f1 R] = deconv(R,D)
```
``` python
print("Drücke den Knopf damit der nächste Rechenschritt angezeigt wird, das Ergebnis wird unter dem Knopf ausgegeben:")
button = widgets.Button(description="Hier drücken :-)")
output = widgets.Output()
%% Cell type:code id: tags:
display(button, output)
``` octave
R = conv(R,[1 0])
[f2 R] = deconv(R,D)
```
i = 0
R = []
fk = []
fo = False
%% Cell type:code id: tags:
def convolve_and_deconvolve(b):
with output:
global i
global R
global fk
global fo
i += 1
``` octave
R = conv(R,[1 0])
[f3 R] = deconv(R,D)
```
#die gleichen Schritte wie in Beispiel 1 (schriftliche Division würde auch hier zum Ergebnis führen):
if i == 1:
f, R = deconvolve(N,D,str(i))
if len(f) == 0:
fo = True
else:
R = convolve(R)
f, R = deconvolve(R,D,str(i))
%% Cell type:code id: tags:
#gebe das momentane Ergebnis für f(k) aus:
fk = list(map(sum, itertools.zip_longest(f, fk, fillvalue=0)))
fk_temp = fk.copy()
if fo == True:
fk_temp.insert(0,0)
print("f(k):"+str(fk_temp))
print("\n-----------------------\n")
``` octave
R = conv(R,[1 0])
[f4 R] = deconv(R,D)
button.on_click(convolve_and_deconvolve)
```
%% Cell type:code id: tags:
%% Output
Drücke den Knopf damit der nächste Rechenschritt angezeigt wird, das Ergebnis wird unter dem Knopf ausgegeben:
``` octave
R = conv(R,[1 0])
[f5 R] = deconv(R,D)
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
- Die Koeffizienten der Folge <code>f</code> entsprechen den Elementen des Zählervektors <code>N</code>.
- Die Koeffizienten der Folge <code>f(k)</code> entsprechen den Elementen des Zählervektors <code>N</code>.
- <code>R</code> wird lediglich = $ 0$.
- Gezeigt wurde ein <span style='color:OrangeRed'>„Finite Impulse Response“ </span>(kurz: <span style='color:OrangeRed'>FIR </span>) Filter.
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #3 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Es gelten weiterhin die Werte aus dem <span style='color:Gray'>Beispiel #1 </span>:
Das Ergebnis (f(k)) aus <span style='color:Gray'>Beispiel #1 </span> wird nun als Funktion geplottet. Es gelten weiterhin die Werte:
$num$=$8z^{-1}$
$den$=$1-3z^{-1}+2z^{-2}$
Wir bestimmen außerdem die Anzahl der Abtastpunkte:
$npoint$ = $20$
Wie zuvor, wird wiederholt die Entfaltung, also Division durchgeführt. Der Code wird in
einer Funktion „ longdiv “ zusammengefasst.
%% Cell type:code id: tags:
``` octave
function y = longdiv(num,den,npoint)
[y(1) R] = deconv(num,den);
for i=2:npoint
num = conv(R,[1 0]);
[f R] = deconv(num,den);
y(i) = f(length(f));
end
end
``` python
def longdiv(N, D, k):
f,R = signal.deconvolve(N,D)
if len(f)==0:
fo = True
else:
fo = False
for i in range(k-1):
R = signal.convolve(R,[1, 0])
f_temp, R = signal.deconvolve(R,D)
f = list(map(sum, itertools.zip_longest(f, f_temp, fillvalue=0)))
if fo == True:
f.insert(0,0)
return(f)
```
%% Cell type:code id: tags:
``` octave
num = [8 0];
den = [1 -3 2];
npoint = 20;
y = longdiv(num, den, npoint);
plot(y)
``` python
num = [8, 0]
den = [1, -3, 2]
npoint = 20
fk = longdiv(num, den, npoint);
print("Werte von f(k) von k=0 bis k=19:"+str(fk))
plt.plot(range(npoint),fk,'blue')
plt.grid()
plt.legend(['f(k)'])
plt.show()
```
%% Output
Werte von f(k) von k=0 bis k=19:[0, 8.0, 24.0, 56.0, 120.0, 248.0, 504.0, 1016.0, 2040.0, 4088.0, 8184.0, 16376.0, 32760.0, 65528.0, 131064.0, 262136.0, 524280.0, 1048568.0, 2097144.0, 4194296.0]
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Das System ist instabil.
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #4 </span>
%% Cell type:code id: tags:
``` octave
pkg load control
``` python
Ts = 0.1
p = -1
npoint = 20
Ts = 0.1;
p = -1;
npoint = 20;
N=[1, 0]
D=[1, -exp(Ts*p)]
y1 = longdiv([1 0],[1 -exp(Ts*p)],npoint);
t =(1:npoint)*Ts-Ts;
yc = longdiv(N, D, npoint);
print("Werte von yc(k) von k=0 bis k="+str(npoint-1)+":"+str(np.around(yc,3)))
plt.plot(np.arange(0.0, npoint*Ts, Ts),yc,'blue')
plt.grid()
plt.legend(['yc'])
plt.show()
```
%% Cell type:code id: tags:
%% Output
``` octave
yc = impulse(tf(1,[1 1]),t);
plot(t,yc,'r',t,y1,'*')
```
Werte von yc(k) von k=0 bis k=19:[1. 0.905 0.819 0.741 0.67 0.607 0.549 0.497 0.449 0.407 0.368 0.333
0.301 0.273 0.247 0.223 0.202 0.183 0.165 0.15 ]
%% Cell type:code id: tags:
``` octave
yh = longdiv(1-exp(Ts*p),[1 -exp(Ts*p)],npoint);
plot(t,yh,'+')
``` python
yh = longdiv(1-exp(Ts*p),[1, -exp(Ts*p)],npoint);
print("Werte von yh(k) von k=0 bis k="+str(npoint-1)+":"+str(np.around(yh,3)))
plt.scatter(np.arange(0.0, npoint*Ts, Ts),yh, marker = '+')
plt.grid()
plt.legend(['yh'])
plt.show()
```
%% Output
Werte von yh(k) von k=0 bis k=19:[0. 0.095 0.086 0.078 0.07 0.064 0.058 0.052 0.047 0.043 0.039 0.035
0.032 0.029 0.026 0.023 0.021 0.019 0.017 0.016]
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #5 </span>
%% Cell type:code id: tags:
``` octave
function [G Ph] = BodeZOH(Ts,ommax,dom)
``` python
def BodeZOH(Ts,om,ommax,dom):
k=0;
kmax = round(ommax/dom);
for k in range(1, kmax):
om[k] = k*dom
Renum = 1 -cos(-om[k]*Ts)
Imnum = -sin(-om[k]*Ts)
Gnum = sqrt(Renum*Renum+Imnum*Imnum)
G[k] = Gnum/(om[k]*Ts)
Ph[k] = atan2(Imnum,Renum)-pi/2
k=0;
kmax = round(ommax/dom);
for k=1:kmax
om(k) = k*dom;
Renum = 1 -cos(-om(k)*Ts);
Imnum = -sin(-om(k)*Ts);
Gnum = sqrt(Renum*Renum+Imnum*Imnum);
G(k) = Gnum/(om(k)*Ts);
Ph(k) = atan2(Imnum,Renum)-pi/2;
end
end
return [G, Ph]
Ts = 0.1;
ommax = 100;
ommax = 200;
dom = 0.1;
om=[0.1:0.1:100];
om=np.arange(1.0, ommax+1, dom)
G=np.arange(1.0, ommax+1, dom)
Ph=np.arange(1.0, ommax+1, dom)
G, Ph = BodeZOH(Ts,om,ommax,dom)
Ph[0]=0.0
plt.plot(om,20*log10(G),'blue')
plt.xscale("log")
plt.grid()
plt.legend(['20*log10(G) (Amplitude)'])
plt.show()
plt.plot(om,Ph*90/pi,'orange')
plt.xscale("log")
plt.grid()
plt.legend(['Ph in Grad (Phase)'])
plt.show()
```
[G Ph] = BodeZOH(Ts,ommax,dom);
%% Output
subplot(2,1,1), semilogx(om,20*log10(G));
subplot(2,1,2), semilogx(om,Ph*90/pi);
```
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Die Amplitude ist im z Bereich kleiner als 20dB.
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #6 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
- Wir stellen drei Verfahren vor, wie in Matlab die Transformation von kontinuierlichen Systemen zu diskreten Systemen erfolgen kann.
- Dabei wird die Matlabfunktion<code>c2d</code> verwendet
- Gegeben sei das System:
$G(s)=\frac{1}{s+1}$ und $T_s=0,1$
%% Cell type:code id: tags:
``` octave
pkg load control
```
``` python
Ts = 0.1
s = control.TransferFunction.s
z = control.TransferFunction.z
print('Abtastzeit: '+str(Ts)+"\n")
%% Cell type:code id: tags:
``` octave
disp('Abtastzeit')
Ts = 0.1
G = (1)/(s + 1)
print('Kontinuierliches System:\n')
print(G)
print("-----------------------------\n")
disp('Kontinuierliche Systeme')
G = tf(1,[1 1])
print('Beispiel 1: ZOH\n')
Gh = cont2dis_to_list([G.num[0][0],G.den[0][0]],Ts,method='zoh')#liefert selbes Ergebnis wie die Matlab-Funktion
Gh = control.TransferFunction(Gh[0],Gh[1],Ts)
print(Gh)
print("-----------------------------\n")
disp('Beispiel 1: ZOH')
Gh = c2d(G,Ts,'zoh')
disp('Beispiel 2: Impulse')
Gz = c2d(G,Ts,'impulse')
print('Beispiel 2: Impulse\n')
Gz = cont2dis_to_list([G.num[0][0],G.den[0][0]],Ts,method='impulse')#!!!!liefert anderes Ergebnis als Matlab-Funktion!!!!Unterschied um den Faktor 10z
Gz = control.TransferFunction(Gz[0],Gz[1],Ts)
print(Gz)
print("-----------------------------\n")
disp('Bespiel 3: Verfahren G(s)/s')
Gi = tf([1],[1 0])
print('Bespiel 3: Verfahren G(s)/s')
print("\nGi:")
Gi = (1)/(s)
print(Gi)
print("Gs = G*Gi:")
Gs = G*Gi
Gsz = c2d(Gs,Ts,'impulse')
Gh = tf([1 -1], [1 0],Ts)
Gf = Gh*Gsz/Ts;
minreal(Gf)
```
print(Gs)
print("Gsz:")
Gsz = cont2dis_to_list([Gs.num[0][0],Gs.den[0][0]],Ts,method='impulse')#!!!!liefert anderes Ergebnis als Matlab-Funktion!!!!Unterschied um den Faktor 10z
Gsz = control.TransferFunction(Gsz[0],Gsz[1],Ts)
print(Gsz)
print("Gh:")
Gh = (z-1)/(z)
Gh.dt = Ts
print(Gh)
print("Gf=Gh*(Gsz/Ts):\n")
Gf = Gh*(Gsz/Ts)
print(Gf)
Gf.minreal()
```
%% Output
Abtastzeit: 0.1
Kontinuierliches System:
1
-----
s + 1
-----------------------------
Beispiel 1: ZOH
0.09516
----------
z - 0.9048
dt = 0.1
-----------------------------
Beispiel 2: Impulse
0.1
----------
z - 0.9048
dt = 0.1
-----------------------------
Bespiel 3: Verfahren G(s)/s
Gi:
1
-
s
Gs = G*Gi:
1
-------
s^2 + s
Gsz:
0.009516
----------------------
z^2 - 1.905 z + 0.9048
dt = 0.1
Gh:
z - 1
-----
z
dt = 0.1
Gf=Gh*(Gsz/Ts):
0.009516 z - 0.009516
--------------------------------
0.1 z^3 - 0.1905 z^2 + 0.09048 z
dt = 0.1
$$\frac{0.09516}{z^2 - 0.9048 z}\quad dt = 0.1$$
TransferFunction(array([0.09516258]), array([ 1. , -0.90483742, 0. ]), 0.1)
%% Cell type:code id: tags:
``` octave
``` python
% Set the Octave Engine to run the simulatio
% Simulation Parameters
pkg load control
addpath("./Octsim");
% Start time
tini = 0;
% End time
tfinal = 5;
% Time Step
dt = 0.001;
% Number of data flows in the schematic
nflows = 5;
% Sampling time for discrete time
Ts = 0.1;
% Instance of the simulation schematic
sc = Schema(tini,tfinal,dt,nflows);
% List of components
c{1} = SinusoidalSignalSource(1,1,2*pi,0);
c{2} = TransferFunction(1,2,1,[1,1]); %G(s)
c{3} = ZOH(1,3,Ts);
c{4} = TransferFunction(3,4,1,[1,1]); %Example 1 - ZOH at the input
c{5} = DTTransferFunction(1,5,[0.09516],[1 -exp(-0.1)],Ts); %Example 2
c{6} = DTTransferFunction(1,6,[0.1 0],[1 -exp(-0.1)],Ts) ;%Example 3
sc.AddListComponents(c);
% Run the schematic and plot
out = sc.Run([ 2 3 4 5 6]);
%plot(out(1,:),out(2,:),out(1,:),out(3,:),out(1,:),out(4,:));
plot(out(1,:),out(2,:),out(1,:),out(4,:),out(1,:),out(5,:), out(1,:),out(6,:));
```
%% Cell type:markdown id: tags:
## <span style='color:Gray'>Beispiel #7 </span>
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Gegeben ist : $G(z)=\frac{z}{z-e^{-\tau}}$, und $U(z)=\frac{z}{z-1}$ mit $\tau=0,1$
Und nehmen wir an: $Y(z)=G(z)U(z)$
%% Cell type:code id: tags:
``` octave
pkg load control
``` python
Ts = 0.1;
#G(z):
Gz=[[1, 0],[1, -exp(-Ts)]]
#%Impuls: longdiv(num,den,100)
y1=longdiv(Gz[0],Gz[1],100);
plt.plot(range(0,100),y1)
plt.grid()
plt.legend(['G(z) => y1(k)'])
plt.show()
```
%% Cell type:code id: tags:
%% Output
``` octave
Ts = 0.1;
Gz = tf([1 0],[1 -exp(-0.1)],Ts);
%Impuls: longdiv(num,den,100)
y1=longdiv(cell2mat(Gz.num),cell2mat(Gz.den),100);
plot(y1);
```
%% Cell type:code id: tags:
``` octave
Uz = tf([1 0],[1 -1],Ts);
Yz = Gz*Uz;
y2=longdiv(cell2mat(Yz.num),cell2mat(Yz.den),100);
plot(y2);
``` python
Uz = [[1, 0],[1, -1]]
Yz = [np.polymul(Uz[0],Gz[0]), np.polymul(Uz[1],Gz[1])]
#Y(z)=G(z)*U(z)
#print(Yz[0])
#print(Yz[1])
y2=longdiv(Yz[0],Yz[1],100);
plt.plot(range(0,100),y2)
plt.grid()
plt.legend(['G(z)*U(z) => y2(k)'])
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
<div style="font-family: 'times'; font-size: 13pt; text-align: justify">
Die Übertragungsfunktion lautet stattdessen:
<br><br>
$G(z)$=$\frac{z(z-e^{-\tau})-(z-1)}{(z-1)(z-e^{-\tau})}$
%% Cell type:code id: tags:
``` octave
Gz2 = tf([1-exp(-0.1) 0],[1 -1-exp(-0.1) exp(-0.1)],Ts);
y3=longdiv(cell2mat(Gz2.num),cell2mat(Gz2.den),100);
plot(y3);
``` python
Gz2 = [[1-exp(-Ts),0],[1,-1-exp(-Ts),exp(-Ts)],Ts]
y3=longdiv(Gz2[0],Gz2[1],100);
plt.plot(range(0,100),y3)
plt.grid()
plt.legend(['y3(k)'])
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` octave
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment