Echantillonnage

1 - Chaine de traitement numérique

Analyse Spectrale au cours de la chaine de traitement :

2 - Définition

$y_e(t) = y(t).\sum_{n=-\infty}^{+\infty} \delta(t-nT_e)$
$y_e(t) = \sum_{n=-\infty}^{+\infty} y(nT_e) \delta(t-nT_e)$

3 - Transformée de Fourier d’un Signal Echantillonné ( Analyse Fréquentielle )

$Y_e(f)= \int_{-\infty}^{+\infty} \color{red}{ y_e(t) } e^{-j2\pi ft}dt$
$Y_e(f)= \int_{-\infty}^{+\infty} \color{red}{ \sum_{n=-\infty}^{+\infty} y(nT_e) \delta(t-nT_e) } e^{-j2\pi ft}dt$

$Y_e(f)= \sum_{n=-\infty}^{+\infty} y(nT_e) \color{green}{ \int_{-\infty}^{+\infty} \delta(t-nT_e) e^{-j2\pi ft}dt }$
$Y_e(f)= \sum_{n=-\infty}^{+\infty} y(nT_e) \color{green}{ e^{-j2\pi nfT_e} }$

Périodicité du Spectre

$Y_e(f+F_e)= \sum_{n=-\infty}^{+\infty} y(nT_e) e^{-j2\pi n(f+F_e)T_e }$
$Y_e(f+F_e)= \sum_{n=-\infty}^{+\infty} y(nT_e) e^{-j2\pi nfT_e } e^{-j2\pi nF_eT_e }$
$Y_e(f+F_e)= \sum_{n=-\infty}^{+\infty} y(nT_e) e^{-j2\pi nfT_e } e^{-j2\pi n }$
Sachant que $e^{-j2\pi n }=1$

$Y_e(f+F_e)= Y_e(f)$

Illustration

Prenons l’exemple de la sinusoïde échantillonnée :

Si j’ajoute progressivement les composantes issues de la série de fourier, nous voyons que la forme mathématique de notre signal échantillonné apparait progressivement en faisant tendre N vers l’infini.
Les composantes résultent de la périodicité du spectre .

Conséquence : Repliement du spectre (aliasing)

Considérons un signal d’entrée un peu bruité.
La périodisation du spectre engendre l’accumulation de ce bruit sur toute la bande passante.
Le signal échantillonné est donc noyé dans le bruit.

Conséquence : Il est nécessaire de filtrer le signal analogique d’entrée.
La fréquence de coupure du filtre doit être Fe/2 ( fréquence de Shannon ).

4 - Transformée de Laplace d’un Signal Echantillonné

$Y_e(s)= \int_{0}^{+\infty} \color{red}{ y_e(t) } e^{-st}dt$
$Y_e(s)= \int_{0}^{+\infty} \color{red}{ \sum_{n=0}^{+\infty} y(nT_e) \delta(t-nT_e) } e^{-st}dt$

$Y_e(s)= \sum_{n=0}^{+\infty} y(nT_e) \color{green}{ \int_{0}^{+\infty} \delta(t-nT_e) e^{-st}dt }$
$Y_e(s)= \sum_{n=0}^{+\infty} y(nT_e) \color{green}{ e^{-s nT_e} }$

Reprenons le cas du signal échelon ; on observe une périodisation de la transformée de Laplace selon la composante $\omega$

5 - Transformée en Z

Reprenons la transformée de Laplace d’un signal échantillonné :

$Y_e(s)= \sum_{n=0}^{+\infty} y(nT_e) e^{-s nT_e}$
posons $z= e^{sT_e}$
Dès lors nous appelons transformée en Z :

$\boxed{ \mathcal{Z}[y_e(t)]=Y(z)= \sum_{n=0}^{+\infty} y(nT_e) z^{-n} }$

Noté également ( sous entendu avec période d’échantillonage $T_e$ :

$\boxed{ Y(z)= \sum_{n=0}^{+\infty} y(n) z^{-n} }$

La transformée en z correspond donc à la transformée de Laplace d’un système échantillonné . Le changement de variable entre s et z nous permettra d’exploiter les calculs de séries mathématiques.

Quant à l’analyse d’un système, elle se fait désormais dans le plan des z :

Comme on est sur une fonction périodique, on retrouve la même chose tous les $\omega _e$, d’où le cercle.

6 - Propriétés de la transformée en Z

6.1 - le retard

soit $y_n=x_{n-k}$ avec $y_n=0$ pour $n<k$
$y_n=x_{n-k}$ est en retard de k échantillons par rapport à y.

$\mathcal{Z}[y_n]=\sum_{n=0}^{+\infty} y_n.z^{-n} = \sum_{n=k}^{+\infty} x_{n-k}.z^{-n}$
$\mathcal{Z}[y_n]=\sum_{m=0}^{+\infty} x_{m}.z^{-m-k} = z^{-k} \sum_{m=0}^{+\infty} x_{m}.z^{-m}$

$\boxed{ \mathcal{Z}[x_{n-k}]=z^{-k}X(z) }$

$z^{-k}$ représente donc un retard de k échantillons.

C’est l’un des points les plus importants : on peut trouver une relation entre les échantillons d’entrée et de sortie ( équation de récurrence ), directement à partir de la transformée en z.

6.2 - Théorème de la valeur finale

$\lim_{z \to 1}(1-z^{-1})Y(z)= \lim_{n \to \infty} yn$

7 - Quelques transformées en Z

• Impulsion : $\mathcal{Z}[\delta(n)]=1$
• Echelon : $\mathcal{Z}[step(n)]=\frac{z}{z-1}=\frac{1}{1-z^{-1}}$
• Exponentiel : $\mathcal{Z}[e^{-\alpha n}]=\frac{z}{z-e^{-\alpha T_e}}=\frac{1}{1-e^{-\alpha T_e}z^{-1}}$
• sinus : $\mathcal{Z}[\sin(n\omega )]=\frac{z. \sin(\omega T_e)}{z^2-2z\cos(\omega T_e)+1}$

8 - Modèle du Bloqueur d’ordre zéro

$b_0(t) = step(t)-step(t-T_e)$
$B_0(s) = \frac{1}{s}-\frac{e^{-sT_e}}{s}$

$B_0(s) = \frac{1-e^{-sT_e}}{s}$

$HB(z)= \mathcal{Z}[B_0(s).H(s)]$
$HB(z)= \mathcal{Z}\left [\frac{1-e^{-sT_e}}{s}.H(s) \right]$
$HB(z)= \mathcal{Z} \left [ 1-e^{-sT_e} \right ].\mathcal{Z} \left [ \frac{H(s)}{s} \right ]$

$\boxed{ HB(z)= (1-z^{-1})Z\left [\frac{H(s)}{s} \right] }$

9 - Approximation numérique de la dérivée

9.1 - Euler

$y(t)=\frac{dx}{dt} = \lim_{\delta t \to 0} \frac{x(t)-x(t-\delta t}{\delta t}$
Pour $T_e$ “petit” :
$y(t)=\frac{dx}{dt} \approx \frac{x(t)-x(t-T_e}{T_e}$
à $t=n.T_e,$
$y(nT_e)=y(n) \approx \frac{x(n)-x(n-1)}{T_e}$

$\boxed{ \frac{dx(t)}{dt} \to \frac{x_n - x_{n-1} }{T_e} }$

En Z :

$Y(z)=\frac{1-z^{-1}}{T_e}.X(z)$ représente donc la dérivée.
On en déduit que s peut être remplacé par $\frac{1-z^{-1}}{T_e}$ pour passer de l’expression analogique en s à l’expression échantillonnée en z.

$\boxed{ s \approx \frac{1-z^{-1}}{T_e} }$

9.2 - Transformation bilinéaire ( méthode des trapèzes )

$x(t)=\frac{dy(t)}{dt}$
$\int_{(n-1)T_e}^{nT_e}x(t)dt \approx T_e\frac{x(n)+x(n-1)}{2}=y(n)-y(n-1)$

$\boxed{ s \approx \frac{2}{T_e} \frac{1-z^{-1}}{1+z^{-1}} }$

10 - Application : Réponse d’un système du 1er ordre

10.1 - Méthode de l’invariant indiciel ( modèle du bloqueur d’ordre 0 )

Système du premier ordre :

$H(s)=\frac{G}{1+ \tau s}$
$HB(z)= (1-z^{-1})\mathcal{Z} \left [\frac{H(s)}{s} \right]$
$HB(z)=(1-z^{-1})\mathcal{Z} \left [ \frac{G}{s(1+ \tau s)} \right ]$

$HB(z)=(1-z^{-1}) \frac{Gz \left (1-e^{-\frac{T_e}{\tau}} \right )}{(z-1) \left (z-e^{-\frac{T_e}{\tau}} \right )}$

$HB(z)= \frac{G \left (1-e^{-\frac{T_e}{\tau}} \right )}{ \left (z-e^{-\frac{T_e}{\tau}} \right )}$

$Y(z)=HB(z).X(z)$

$Y(z)= \frac{G \left (1-e^{-\frac{T_e}{\tau}} \right )}{ \left (z-e^{-\frac{T_e}{\tau}} \right )} .\frac{z}{z-1}$

Transformation inverse (utilisation des tables) :

$\boxed{ y_n = G.(1-e^{-n.\frac{T_e}{\tau}}) }$

ou détermination de l’équation de récurrence :

$\frac{Y(z)}{X(z)}= \frac{G \left (1-e^{-\frac{T_e}{\tau}} \right )}{ \left (z-e^{-\frac{T_e}{\tau}} \right )}$

Passage en $z^{-1}$ :

$\frac{Y(z)}{X(z)}= \frac{G \left (1-e^{-\frac{T_e}{\tau}} \right )}{ \left (1-z^{-1}e^{-\frac{T_e}{\tau}} \right )}$

$\left (1-z^{-1}e^{-\frac{T_e}{\tau}} \right ).Y(z)= G \left (1-e^{-\frac{T_e}{\tau}} \right ).X(z)$
$Y(z)-z^{-1}Y(z)e^{-\frac{T_e}{\tau}} = G \left (1-e^{-\frac{T_e}{\tau}} \right ).X(z)$

$Y(z) = z^{-1}Y(z)e^{-\frac{T_e}{\tau}} + G \left (1-e^{-\frac{T_e}{\tau}} \right ).X(z)$

$\boxed{ y_n=e^{-\frac{T_e}{\tau}} y_{n-1} + G \left (1-e^{-\frac{T_e}{\tau}} \right ) x_n }$

$x_n = 1$ pour $n \geqslant 0$
$y_0 = 0$

10.2 - Approximation de la dérivée par la méthode d’Euler

$\tau \frac{dy(t)}{dt}+y(t)=G.x(t)$

Dans le cas d’un système échantillonné :

$y \to y_n$
$x \to x_n$
$\frac{dy(t)}{dt} \to \frac{y_n - y_{n-1} }{T_e}$ ( approximation méthode d’Euler )

L’équation devient :

$\tau \frac{y_n - y_{n-1} }{T_e}+y_n=G.x_n$

On en déduit l’équation de récurrence :

$\boxed{ y_n=\frac{\tau}{\tau + T_e}.y_{n-1} + \frac{G T_e}{\tau+T_e}.x_n }$

$x_n = 1$ pour $n \geqslant 0$
$y_0 = 0$

Cette équation de récurrence peut être déduite à partir de la fonction de transfert :

$H(s) = \frac{G}{1+ \tau s}$
$s \approx \frac{1-z^{-1}}{T_e}$

$H(z) = \frac{Y(z)}{X(z)} = \frac{G}{1+ \tau \frac{1-z^{-1}}{T_e} }$

$H(z) = \frac{G}{1+ \tau \frac{1-z^{-1}}{T_e} }$

$Y(z)=\frac{\tau}{\tau + T_e}.z_{-1}.Y(z) + \frac{G T_e}{\tau+T_e}.X(z)$

On retrouve :

$y_n=\frac{\tau}{\tau + T_e}.y_{n-1} + \frac{G T_e}{\tau+T_e}.x_n$

REMARQUE : On peut trouver l’expression de y en passant par Z :

Repartons de l’équation différentielle “numérisée” :

$y_n - K_2.y_{n-1}=G.K_1.x_n$ avec $K_1=\frac{T_e}{\tau+T_e}$ et $K_2=1-K_1$

$\mathcal{Z}\left [ y_n - K_2.y_{n-1} \right ] =\mathcal{Z} \left [ G.K_1.x_n \right ]$

$Y(z)-K_2z^{-1}Y(z)=G.K_1.X(z)$

$\boxed{ \frac{Y(z)}{X(z)}=\frac{G.K_1}{1-K_2z^{-1}} }$

Réponse Indicielle ( echelon)

$X(z)=\frac{1}{1-z^{-1}}$
$Y(z)=\frac{G.K}{1-K_2z^{-1}}.\frac{1}{1-z^{-1}}$

$\boxed{ Y(z)=GK_1\frac{z^2}{[z-K_2][z-1]} }$

Décomposition en éléments simples :

$\frac{Y(z)}{z}=\frac{A}{z-K_2}+\frac{B}{z-1}$

$A = \lim_{z \to K_2} (z-K_2)\frac{Y(z)}{z}=-G.K_2$
$B = \lim_{z \to 1} (z-1)\frac{Y(z)}{z}=G$

$Y(z) = \frac{Gz}{z-1}+\frac{-G.K_2.z}{z-K_2}$

Utilisation des tables pour trouver $y_n$

$\boxed{ y_n = G(1-K_2^{n}) }$

10.3 - Approximation de la dérivée par la transformation bilinéaire

$s \approx \frac{2}{T_e} \frac{1-z^{-1}}{1+z^{-1}}$

$H(s) = \frac{G}{1+ \tau s}$

$H(z) = \frac{G}{1+\tau \left [ \frac{2}{T_e} \frac{1-z^{-1}}{1+z^{-1}} \right ]}$

$H(z) = \frac{G(1+z)}{K_1.z - K_2 }$
avec $K_1=1+\frac{2\tau}{T_e}$ et $K_2=\frac{2\tau}{T_e} - 1$

$H(z) = \frac{Y(z)}{X(z)} = \frac{G}{K_1} . \frac{1+z^{-1}}{1 - K_2.z^{-1} }$

Equation de récurrence :

$\boxed{ y_n=\frac{K_2}{K_1}.y_{n-1} + \frac{G}{K_1}x_n + \frac{G}{K_1}x_{n-1} }$

10.4 - Résumé

11 - Script Scilab / octave-matlab pour générer les figures ci dessus :

figure(1)


//==============================================================================
// CHAINE DE TRAITEMENT NUMERIQUE
//==============================================================================

Npp=40 ; a=1; Nd=80; N=Npp*Nd;
Fe=16e3; Te=1/Fe; Fc=Fe*Npp; Tc=1/Fc; F=1000;
t=Tc*(0:N-1);
ve=5*sin(2*%pi*F*t);
d=[ones(1, a), zeros(1, (Npp-a))];
imp=matrix(d' *ones(1, Nd ), 1 , N);
se=ve .* imp;
// réponse impulsionnelle du bloqueur
b=ones(1, Npp);
// signal en sortie du bloqueur, réajusté sur N points
v1= convol(b,se);
v1=v1(1 : N);
// calcul du filtre passe bas, butteworth d'ordre 8
G=analpf(8,'butt',[0 0],2*%pi*(0.4*Fe));
f=Fc/N*[0 : N-1];
GAIN=freq(G(2),G(3),%i*2*%pi*f);
G=[GAIN(1:N/2), conj(GAIN(N/2+1:-1:2))];
hpb=fft(G, 1);
// filtrage
v2=convol(hpb, v1); v2=v2(1:N);


//==============================================================================
// Display
//==============================================================================

t_min=0
t_max=0.004
y_min=-6
y_max=6

subplot(4,2,1)
plot2d(t,ve); 
    p = get("hdl")
    p.children.thickness = 3;
    p.children.foreground = 5; 
    a=gca();
    a.data_bounds=[t_min,t_max,y_min,y_max];
    a.x_location = "origin";
    a.y_location = "origin"; 
    xtitle("temp signal","t (s)","amp (V)");
plot2d2(t,imp);
xtitle("Input signal, Dirac","t (s)","amplitude (V)")


subplot(4,2,2)
Fs=640000
S1=1/N*fft(ve,-1);
S=[S1(N/2+1:N), S1(1:N/2)];
f=(Fs/N)*[-N/2 : N/2-1];
plot2d3(f,abs(S),5)
    p = get("hdl");
    p.font_foreground=2;
    p.font_style=5;
    p.children.thickness = 3;
    p.children.foreground = 3; 
    a=gca();
    a.x_location = "origin";
    a.y_location = "origin";
    a.data_bounds=[-3*F,3*F,0,2.5]; // f_min,f_max,yf_min,yf_max
    xtitle("spectrum analysis"," fréquency (Hz)","amp (Vs)");



//------------------------------------------------------------------------------
subplot(4,2,3)
plot2d2(t,se);
    p = get("hdl")
    p.children.thickness = 3;
    p.children.foreground = 9; 
    a=gca();
    a.data_bounds=[t_min,t_max,y_min,y_max];
    a.x_location = "origin";
    a.y_location = "origin"; 
    xtitle("temp signal","t (s)","amp (V)");
xtitle("Sampled signal","t (s)","amplitude (V)");


subplot(4,2,4)
Fs=640000
S1=1/N*fft(se,-1);
S=[S1(N/2+1:N), S1(1:N/2)];
f=(Fs/N)*[-N/2 : N/2-1];
plot2d3(f,abs(S),5)
    p = get("hdl");
    p.font_foreground=2;
    p.font_style=5;
    p.children.thickness = 3;
    p.children.foreground = 3; 
    a=gca();
    a.x_location = "origin";
    a.y_location = "origin";
    a.data_bounds=[-2*Fe,2*Fe,0,0.1]; // f_min,f_max,yf_min,yf_max
    xtitle("spectrum analysis"," fréquency (Hz)","amp (Vs)");

//------------------------------------------------------------------------------
subplot(4,2,5)
plot2d2(t,v1);
    p = get("hdl")
    p.children.thickness = 3;
    p.children.foreground = 5; 
    a=gca();
    a.data_bounds=[t_min,t_max,y_min,y_max];
    a.x_location = "origin";
    a.y_location = "origin";
xtitle("Bloqued Signal","t (s)","amplitude (V)")
 


subplot(4,2,6)
Fs=640000
S1=1/N*fft(v1,-1);
S=[S1(N/2+1:N), S1(1:N/2)];
f=(Fs/N)*[-N/2 : N/2-1];
plot2d3(f,abs(S),5)
    p = get("hdl");
    p.font_foreground=2;
    p.font_style=5;
    p.children.thickness = 3;
    p.children.foreground = 3; 
    a=gca();
    a.x_location = "origin";
    a.y_location = "origin";
    a.data_bounds=[-10*Fe,10*Fe,0,1]; // f_min,f_max,yf_min,yf_max
    xtitle("spectrum analysis"," fréquency (Hz)","amp (Vs)");


//------------------------------------------------------------------------------
subplot(4,2,7)
plot2d2(t,v2);
    p = get("hdl")
    p.children.thickness = 3;
    p.children.foreground = 5; 
    a=gca();
    a.data_bounds=[t_min,t_max,y_min,y_max];
    a.x_location = "origin";
    a.y_location = "origin"; 
xtitle("Analog Filtered output signal","t (s)","amplitude (V)");



subplot(4,2,8)
Fs=640000
S1=1/N*fft(v2,-1);
S=[S1(N/2+1:N), S1(1:N/2)];
f=(Fs/N)*[-N/2 : N/2-1];
plot2d3(f,abs(S),5)
    p = get("hdl");
    p.font_foreground=2;
    p.font_style=5;
    p.children.thickness = 3;
    p.children.foreground = 3; 
    a=gca();
    a.x_location = "origin";
    a.y_location = "origin";
    a.data_bounds=[-2*Fe,2*Fe,0,2.5]; // f_min,f_max,yf_min,yf_max
    xtitle("spectrum analysis"," fréquency (Hz)","amp (Vs)");


//==============================================================================
//              PERIODISATION DU SPECTRE : ILLUSTRATION
//==============================================================================
figure(2)
t=0:0.0001:3;
F=1
Fe=10
omega=2*%pi*F
N_vect=[3 10 50 1000]

for i=1:4
//-----------------------------------------------------------------------------
        subplot(2,2,i)
        N=N_vect(i)
        s = 0//cos(omega*t)
        for k=0:N
               
                s= s + (1/N)*(0.5*cos(2*%pi*(k*Fe-F)*t)+0.5*cos(2*%pi*(k*Fe+F)*t));
        end
        plot2d(t,s)
        p = get("hdl");
        p.children.thickness = 2;
        p.children.foreground = 5; 
        title([['$y(t)=cos( \omega nT_e)$'], ' deduced from Fourier Serie with N=',string(N),' Terms'],'FontSize',4);
        xlabel('t','FontSize',4);  
        ylabel('y(t)','FontSize',4);  

//-----------------------------------------------------------------------------
end


//==============================================================================
//                  SIGNAL CARRE ECHANTILLONNE
//==============================================================================
N = 1000;       
nd=100
Te=10e-3 
Fe=1/Te
t=Te*(-N/2: N/2-1);

f=1
T=1/f
tau=0.5

t_min=-4
t_max=4
y_min=-0.1
y_max=1.2

f_min=-50
f_max=50
yf_min=0
yf_max=0.7
//------------------------------------------------------------------------------
subplot(2,2,1)
duty=100*(tau/T)
square =0.5 *( squarewave(2*%pi*f*t,duty) + 1)
plot(t,square)
    p = get("hdl");
    p.children.thickness = 4;
    p.children.foreground = 5; 
    a=gca();
    a.data_bounds=[t_min,t_max,y_min,y_max];
    a.x_location = "origin";
    a.y_location = "origin"; 
    xtitle("temp signal","t (s)","amp (V)");

subplot(2,2,2)
S1=1/N*fft(square,-1);
S=[S1(N/2+1:N), S1(1:N/2)];
f=(Fe/N)*[-N/2 : N/2-1];
plot2d2(f,abs(S))
    p = get("hdl");
    p.font_foreground=2;
    p.font_style=5;
    p.children.thickness = 3;
    p.children.foreground = 3; 
    a=gca();
    a.x_location = "origin";
    a.y_location = "origin";
    a.data_bounds=[f_min,f_max,yf_min,yf_max];
    xtitle("spectrum analysis"," fréquency (Hz)","amp (Vs)");

//------------------------------------------------------------------------------
subplot(2,2,3)
    Nd=100
    Npp=N/Nd
    d=[ones(1, a), zeros(1, (Npp-a))];
    imp=matrix(d' *ones(1, Nd ), 1 , N);
    se=square .* imp;
    plot2d2(t,se);
    p = get("hdl")
    p.children.thickness = 3;
    p.children.foreground = 9; 
    a=gca();
    a.data_bounds=[t_min,t_max,y_min,y_max];
    a.x_location = "origin";
    a.y_location = "origin"; 
    xtitle("Sampled signal","t (s)","amplitude (V)");

subplot(2,2,4)
S1=10*1/N*fft(se,-1);
S=[S1(N/2+1:N), S1(1:N/2)];
f=(Fe/N)*[-N/2 : N/2-1];
plot2d2(f,abs(S))
    p = get("hdl");
    p.font_foreground=2;
    p.font_style=5;
    p.children.thickness = 3;
    p.children.foreground = 3; 
    a=gca();
    a.x_location = "origin";
    a.y_location = "origin";
    a.data_bounds=[f_min,f_max,yf_min,yf_max];
    xtitle("spectrum analysis"," fréquency (Hz)","amp (Vs)");

G=10
tau=1
Te=0.1
N=50

v_time=0:Te:(N-1)*Te;

% Analog Output
plot(v_time,G*(1-exp(-v_time/tau)));
hold on;

%-------------------------------------------------------------------------------
%         FIRST METHOD : Zero Order Blocker
%-------------------------------------------------------------------------------
% Output
for n=0:N-1
y(n+1)=G*(1-exp(-n*Te));
end 

stem(v_time,y)
hold on

% recurrence relation
yn_1=0
yn(1)=0
for n=1:N-1
yn(n+1)=yn_1*exp(-Te/tau)+G*(1-exp(-Te/tau));
yn_1=yn(n+1);
end

# stem(v_time,yn)
# hold on

%-------------------------------------------------------------------------------
%         SECOND METHOD : Euler Approx
%-------------------------------------------------------------------------------

K1=Te/(tau+Te)
K2=1-K1

for n=0:N-1
y2(n+1)=G*(1-K2^(n));
end 

stem(v_time,y2)
hold on

%-------------------------------------------------------------------------------
%         THIRD METHOD : Bilinear Approx
%-------------------------------------------------------------------------------
K1=1+2*tau/Te
K2=2*tau/Te-1

% recurrence relation
yn_1=0
ybn(1)=0
for n=1:N-1
ybn(n+1)=yn_1*(K2/K1)+G/K1+G/K1;
yn_1=ybn(n+1);
end

stem(v_time,ybn)
hold on


%-------------------------------------------------------------------------------

legend('analog','zero order blocker','euler approx','bilinear transform')


# print -dsvg ../illustrations/temp_fo.svg