Circuiti lineari frazionali

Relatore:

Massimo Panella

Controrelatore:

Frank Silvio Marzano

Roberto Patrizi

Studio di

Circuiti Lineari Frazionali nel dominio tempo discreto

Roberto Patrizi

Stampato a Frosinone, 18 luglio 2009

e-mail: [email protected]

Tesi specialistica del corso di:Circuiti Intelligenti, Prof. Massimo Panella.Giovedì 23 Luglio 2009, Sapienza, Facoltà di IngegneriaVia Eudossiana 18 Roma

Composto in LATEX

[email protected]

I N D I C E

PREFAZIONE iii

RINGRAZIAMENTI v

TERMINI E SIMBOLI vii

ACRONIMI ix

1 B A S I M AT E M AT I C H E 1

1.1 FUNZIONE GAMMA 1

1.2 DIFFERINTEGRALI 3

1.2.1 Percorso storico 3

1.2.2 Differintegrali di Grünwald-Letnikov 4

1.2.3 Differintegrali di Riemann-Liouville 6

1.2.4 Differintegrali di Courant e Hilbert 7

1.2.5 Differintegrali di funzioni elementari 7

1.3 PROPRIETÀ DIFFERINTEGRALI 9

1.3.1 Differintegrali multipli 9

1.3.2 Differintegrali del prodotto di funzioni 11

1.3.3 Differintegrali di funzioni composte 11

1.4 TRASFORMATA DI LAPLACE 12

1.4.1 Convoluzione e prodotto 13

1.4.2 Proprietà della derivata 13

1.4.3 Teorema del valore iniziale e finale 15

1.5 SUPERFICI DI RIEMANN 15

1.5.1 Radice 15

1.5.2 Logaritmo 16

1.5.3 Rappresentazione grafica 17

1.6 FUNZIONE E DI MITTAG-LEFFLER 18

1.6.1 Trasformata di Laplace di E 18

2 S I S T E M I A N A LO G I C I 21

2.1 DEFINIZIONI GENERALI 21

2.2 SISTEMI ELEMENTARI 22

2.2.1 Termine monomio 22

2.2.2 Termine binomio 24

2.2.3 Funzioni più complesse 25

2.3 STABILITÀ 26

2.4 APPROSSIMAZIONI 27

2.4.1 Espansione in frazioni continue 28

2.4.2 Metodo di Carlson 32

2.4.3 Metodo di Chareff 33

2.4.4 Metodo di Oustaloup 35

2.5 SISTEMI DI CONTROLLO 35

2.5.1 Controllore lead-lag 36

2.5.2 Controllori di tipo PID 37

2.5.3 Controllori di tipo TID frazionari 37

2.5.4 Controllori CRONE 38

i

ii Indice

3 S I S T E M I T E M P O - D I S C R E T I 39

3.1 INTRODUZIONE 39

3.2 RITARDO ELEMENTARE ZA40

3.2.1 Differenze frazionarie 43

3.2.2 Risposta in frequenza 45

3.2.3 Trasformata Z 45


3.3.1 Modelli MA 47

3.3.2 Modello completo 47

3.3.3 Inversione fratti semplici 48

3.4 APPROCCIO ALTERNATIVO 49

3.4.1 Discretizzazione delle derivate non intere 50

3.4.2 Sistemi descritti da equazioni alle differenze finite 51

3.4.3 Algoritmi per equazioni alle differenze 52

3.5 APPROSSIMAZIONI DISCRETE 54

3.5.1 Necessità delle approssimazioni 54

3.5.2 Metodo dell’errore quadratico minimo 56

3.5.3 Metodo di Eulero-Grünwald 57

3.5.4 Metodo di Tustin 59

3.5.5 Metodo di Al-Alaoui 62

3.5.6 Padé, Prony e Shanks 62

4 A P P L I C A Z I O N I E T E C N O LO G I E 65

4.1 ELEMENTI FRAZIONARI CONCENTRATI 65

4.1.1 Condensatore frazionario 65

4.1.2 Induttore frazionario 67

4.2 LINEE DI TRASMISSIONE 68

4.3 PRESTAZIONI ALGORITMI NUMERICI 70

4.3.1 Soluzione analitica 71

4.3.2 Metodi alle differenze 72

4.3.3 Riepilogo risultati 76

4.4 ESEMPI APPLICATIVI 78

4.4.1 Oscillatore smorzato 78

4.4.2 Codifica del parlato 80

4.4.3 Conclusioni 81

A TO O L B O X N I N T E G E R 83

A.1 PROGRAMMAZIONE AD OGGETTI 83

A.2 NINTEGER: COMANDI COMUNI 84

BIBLIOGRAFIA 92

INDICE ANALITICO 94

P R E FA Z I O N E

“This is an apparent paradox from wich,one day, useful consequences will be drawn”

Leibniz, 30 settembre 16951

Di cosa parla questo testo?

In questo scritto saranno descritte e studiate le basi teoriche all’originedi una nuova classe di sistemi elettronici, cioè i sistemi lineari di ordinenon intero. Un sistema di ordine non intero è descritto da potenze nonintere nel dominio di Laplace, come ad esempio il semplice sistemadescritto dalla funzione di trasferimento sα, con α ∈ R. Le equazionicorrispondenti nel dominio del tempo assumono la forma di derivatedi ordine non intero, come mostrato nel § 1.4.2 nella pagina 13.

Perché si studiano questi sistemi, quali sono i vantaggi?

Da qualche anno è cresciuto l’interesse verso sistemi elettronici di que-sto tipo, grazie ai risultati incoraggianti, come dimostrano le prestazio-ni, nettamente superiori, di controllori realizzati con le nuove teorie,rispetto ai controllori standard. Molteplici articoli guardano a questisistemi come ai sistemi del futuro: Ortigueira [30] sottotitola il suo ar-ticolo, che dà una visione d’insieme dello studio dei sistemi frazionari,con le parole:

The 21st Century systems;

Secondo Chen [12] i sistemi frazionari saranno universalmente dif-fusi, per un semplice motivo, che detto con le sue parole:

Using several real world examples, we further argue that,fractional order control is ubiquitous when the dynamicsystem is of distributed parameter nature.

Quindi, ciascun sistema a costanti distribuite, è un sistema frazionarioper sua natura, ed è quindi ben descritto dalle teorie trattate in questotesto.

A che punto sono gli studi?

Le basi dei sistemi di ordine non intero risalgono a più di 300 anni fa,quando stavano venendo alla luce le teorie dell’analisi infinitesimale.Se ha senso calcolare le derivate di ordine intero di una funzione, chesignificato potrebbe avere una derivata di ordine non intero? A partireda questa domanda posta da Leibniz a l’Hôpital, hanno avuto originegli studi sulle teorie matematiche alla base dei sistemi che vedremo inquesto testo, che ormai, con circa 300 anni alle spalle, sono sufficien-temente sviluppate per l’utilizzo. Tuttavia, è solo in tempi recenti chequesti studi hanno cominciato a sfociare in interessanti applicazioni

1 Tratto dalla lettera di risposta del 30 settembre 1695 di Leibniz a l’Hôpital, riguardo alsignificato della scrittura dn f /dxn, se ad esempio si avesse n = 1/2. Disponibile su [3,pagina vii], anche sul sito Google libri.

iii

http://books.google.it/books?id=LhkO83ZioQkC

iv PREFAZIONE

pratiche, dopo essere state a lungo considerate come un esercizio dipensiero.

L’ambito dei controlli è stato il primo ad essere affrontato, attual-mente è il più esteso, sebbene ancora non siano stati sviluppati tuttigli strumenti frazionari di analisi e progettazione, disponibili inveceper i sistemi di ordine intero. Per quanto riguarda i sistemi digitaligli studi sono più carenti, spesso limitati ad implementazioni digitalidi sistemi di controllo. Tuttavia i sistemi digitali di ordine non intero,si comportano come sistemi a ritardo frazionario, appoggiandosi allateoria nota del ricampionamento digitale dei sistemi.

Come documentarsi sull’argomento? Perché questo testo?

Attualmente, per quanto ci è noto, non esistono libri esaustivi riguar-danti la teoria e le applicazioni ai sistemi elettronici di ordine nonintero, a parte forse il recente [16], consultato poco prima di ultimarequesto lavoro. La maggior parte degli studi è disponibile sotto formadi articoli, che spesso trattano un problema specifico.

In questo testo si è cercato di dare una panoramica dello stato attua-le della ricerca, illustrando tutte le conoscenze di base necessarie perintraprendere uno studio dei sistemi razionali non interi. La teoria èstata ricostruita a partire proprio dagli articoli reperibili, pressoché uni-ca e molteplice fonte di informazioni, integrando diverse fonti, esplici-tando calcoli matematici, affiancando grafici e script Matlab, e spessoanche aggiungendo considerazioni personali. Si è quindi cercato diorganizzare in maniera coerente e sistematica, in un unico testo, tuttele informazioni basilari per lo studio di tali sistemi, tenendo presenteun target di lettori che, come l’autore, possiedono una base culturalescientifica data dai corsi di ingegneria.

Come è strutturato questo testo?

Nelle pagine seguenti, si raccolgono le nozioni fondamentali per faci-litare la lettura di questo testo: in particolare ci stiamo riferendo allaterminologia relativa ai nuovi concetti esposti, e alla notazione usata,disponibili nella prossima sezione, dal titolo Termini e simboli; nellasezione successiva, acronimi, sono raccolte tutte le abbreviazioni del-le sigle più comuni, frequentemente utilizzate sia in questo testo chenella maggior parte degli articoli consultati.

Lo studio vero e proprio inizia dal primo capitolo, nel quale si for-niscono tutti gli strumenti matematici alla base dei nuovi sistemi, inparticolare le derivate frazionarie e la relativa teoria delle trasforma-te di ordine non intero, insieme ad altri concetti ausiliari ma non perquesto meno importanti, come la teoria delle superfici di Riemann e lafunzione gamma di Eulero.

Nel secondo capitolo, saranno affrontati i sistemi nel dominio con-tinuo del tempo o di Laplace. Lo scopo è quello di capire come siaffronta lo studio dei sistemi non interi tempo-continui, ricalcandoquindi tutte le procedure note per i sistemi di ordine intero. Le dif-ferenze che emergono, comportano la necessità di introdurre nuovimetodi, in particolare si rendono necessari metodi di approssimazioneper la realizzazione dei sistemi proposti. Uno sguardo infine ai sistemidi controllo, chiude il capitolo.

PREFAZIONE v

Nel terzo capitolo saranno affrontati i sistemi nel dominio discre-to. La discretizzazione dei sistemi analogici pone ulteriori problemi,l’assenza della trasformata Z inversa di un termine del tipo zα inparticolare, richiede una nuova formulazione di antitrasformata Z !

È in questo ambito che si inserisce lo studio specifico di questo testo,che propone un modo alternativo (§ 3.4 nella pagina 49) di trattare i si-stemi digitali, basato sostanzialmente sulla più complessa definizionedi derivata frazionaria di Grünwald-Letnikov (G-L) piuttosto che sulladefinizione, per alcuni aspetti più pratica, data da Riemann e Liouville.È Questa la sezione caratterizzante del lavoro, quello “Studio di circuitilineari frazionali nel dominio tempo discreto” descritto dal titolo. La se-zione sviluppata in modo autonomo, è stata poi ritrovata poco primache questo testo fosse completato, nel recente [16], come metodo dicalcolo.

Nel quarto capitolo invece cercheremo di capire il modo per realiz-zare2 concretamente i sistemi proposti. Saranno analizzati i due princi-pali approcci esistenti, l’approssimazione dei sistemi di ordine non in-tero con componenti discreti di ordine intero, oppure la realizzazionead-hoc di componenti elettronici di ordine non intero.

Per concludere saranno considerati alcuni esempi applicativi, in par-ticolare l’oscillatore di ordine non intero ci permette di applicare su uncaso concreto la teoria sviluppata per i sistemi analogici, mentre peri sistemi di sintesi del parlato, molto più complessi ma di maggioreinteresse pratico e concreto, cercheremo di capire in che modo è pos-sibile migliorare le prestazioni degli algoritmi esistenti focalizzandociproprio sulla sintesi delle vocali!

Parallelamente agli studi teorici, si stanno sviluppando dei pacchettiapplicativi che permettano di studiare i sitemi di ordine non intero. Inparticolare per Matlab, le lacune del controlsystem toolbox sui site-mi di ordine non intero, stanno venendo gradualmente colmate dallacomparsa di pacchetti specifici, dedicati ai nuovi sistemi. Una panora-mica è disponibile nell’appendice A, che fornisce indicazioni utili perapplicare da subito le teorie presentate, ed una serie di riferimenti gra-zie ai quali è possibile trovare rapidamente vari script già pronti perMatlab.

R I N G R A Z I A M E N T I

Non solo un lavoro come questo documento, ma tutta la carrierauniversitaria non sarebbe mai stata possibile senza un valido suppor-to, sia pratico che, soprattutto, morale. Pertanto desidero ringrazia-re in primo luogo i professori della facoltà di ingegneria, che hannosempre mostrato disponibilità e flessibilità, rendendo in questo modo,più semplice il percorso universitario di noi studenti, ed in particolarmodo vorrei ringraziare per questo il mio relatore, il prof. Panella.

Un secondo ringraziamento è per i miei genitori, che mi hanno per-messo di studiare e mi hanno sempre sostenuto ed incoraggiato. Rin-grazio mio fratello Domenico che mi ha spesso aiutato negli studi, così

2 In questo caso stiamo parlando della realizzazione dei sistemi analogici, dal momentoche i sistemi discreti sono realizzabili dall’algoritmo che li definisce, come mostrato nel4° capitolo sui sistemi digitali

vi PREFAZIONE

come i miei amici dell’università che non mi hanno mai negato unaiuto, in particolare Davide, Giovanni e Stefano.

Ringrazio Andrea, che nei viaggi in treno tra Roma e Frosinone haspesso chiarito i miei dubbi e mi ha ascoltato, e gli altri amici di viaggioche hanno pazientemente sopportato la mia voce; mia cugina Lorettache mi ha prestato la sua connessione per scaricare decine di artico-li ed email, e ringrazio tutti quelli che, per brevità, non ho potutoringraziare esplicitamente, compreso me stesso.

Infine, un ringraziamento speciale a Sara, che è stata la mia vitaed il mio mondo d’evasione in questi anni che altrimenti sarebberostati troppo incentrati sul solo studio, che mi ha regalato il sorriso,l’allegria, la vitalità le emozioni e molto di più, in tantissime occasioniche grazie a lei si sono trasformate in momenti speciali e preziosi dacustodire e conservare, e soprattutto che è stata in grado di scopriree farmi conoscere, come nessuno, il lato più umano e remoto del mioanimo! Grazie davvero, buona lettura!

T E R M I N I . . .

Per facilitare la lettura di questo testo e degli articoli riguardanti gliargomenti trattati, si è pensato di raccogliere gli acronimi più diffusiin una comoda tabella, disponibile alle pagine seguenti. Si riportanoinoltre, qui di seguito, i termini relativi ai nuovi sviluppi tecnologicie teorici, ed un elenco di convenzioni matematiche adottate in questotesto.

Il primo e forse il più abusato termine da introdurre è il termine diffe-rintegrale. Quando si parla di derivate ed integrali, si considerano dueconcetti distinti. Nel caso non intero tale differenza scompare! Infatti,in tal caso, considerando ad esempio una derivata, un qualsiasi ordinedi derivazione α ∈ sarà ammissibile, anche negativo! Come è ovvio im-maginare, un ordine di derivazione negativo dà luogo ad un integrale,per cui la differenza tra i due operatore svanisce di fatto. Proprio perquesto viene coniato il il termine differintegrale, che generalizza sia ilconcetto di derivata che quello di integrale. In letteratura, il termineinglese corrispondente è differintegral. La notazione dei differintegraliè varia tanto quanto quella delle derivate.

Si osserva poi, che per brevità, con il nome di sistemi frazionari, sa-ranno indicati in generale i sistemi di ordine non intero, sebbene iltermine possa far pensare ad una minore generalità.

Resta da chiarire cosa si intende per ordine di un sistema: nellarappresentazione nel dominio di Laplace, l’ordine è l’esponente del-la variabile di Laplace comunemente indicata con s; ciò equivale, neldominio del tempo, ad un’equazione in cui compaiono derivate dellostesso ordine dell’esponente della s.

L’impedenza, che comunemente indica la resistenza al passaggio dicorrente di un elemento, come valore complesso dalla variabile di La-place s, prenderà il nome di frattanza quando l’esponente della s èfrazionario. L’equivalente corrispettivo inglese che più spesso ricorrein letteratura è fractance.

Per le impedenze della forma 1/sα, sappiamo che, nel caso ordina-rio con α = 1, stiamo trattando con le impedenze di un condensato-re o capacitore; non c’è un termine altrettanto chiaro nel caso in cuiα ∈ (0, 1). In letteratura si trova spessissimo il termine fractor, che initaliano potrebbe essere reso con frattore, ma che più spesso in questotesto sarà indicato con capacità frazionaria, analogamente per l’induttorefrazionario.

. . . E S I M B O L I

I differintegrali presentano una grande varietà di possibili scritture,un po’ come le derivate ordinarie, alle quali si aggiungono le conven-zioni diffuse per la notazione di ordine non intero. È bene chiarire cheun differintegrale è definito dall’ordine, dall’estremo inferiore considera-to e dall’estremo superiore che coincide con la variabile della funzionedifferintegrata.

vii

viii TERMINI E SIMBOLI

Tabella 1: Elenco delle convenzioni matematiche e dei simboli adottati inquesto testo

Notazione Descrizione

aDαx f Simbolo compatto operatore differintegrale

dα f[d(x− a)]α

Simbolo analitico operatore differintegrale

f (α) Funzione differintegrata

s Variabile complessa di Laplace

j Unità immaginaria

Z {} , L {} Trasformata Zeta e di Laplace. Conesponente −1 si hanno le trasformazioniinverse

N, Z, Q, R, C Insiemi numerici.

⊗ Prodotto di convoluzione

a Vettore. I caratteri in neretto indicano quan-tità vettoriali. Se accostati ad un numerotra parentesi quadre, il numero (od equiva-lentemente il parametro al suo posto) indi-ca l’elemento del vettore considerato, cioè sea = {1, 4, 2, 6}, a[3] = 2.

L’ORDINE tipicamente indicato con n se intero, con α altrimenti, è no-to dal caso intero. L’unica osservazione è che α potrebbe essereun elemento dell’insieme dei numeri complessi C, reali R, o ra-zionali Q, in questo ultimo caso potrebbe anche essere indicatocome 1/q se α ∈ (0, 1).L’ESTREMO INFERIORE a è un concetto noto nel caso del calcolo inte-grale, ma non è tipicamente presente se si considerano le deri-vate. Se però ricordiamo che un differintegrale comprende en-trambi gli operatori, di derivata ed integrale, si osserva comel’estremo inferiore sia necessario sempre, e non pone probleminel calcolo della derivata che se è definita a partire da −∞, saràdefinita a partire da a.L’ESTREMO SUPERIORE x è la variabile rispetto alla quale stiamo cal-colando la derivata, oppure è l’estremo superiore di integrazione,se stiamo calcolando l’integrale rispetto ad una variabile ausilia-ria.

Le notazioni possibili sono elencate nelle prime righe della tabella 1,in cui f è una qualunque funzione differintegrabile. Si osservi che spes-so nella definizione analitica dell’operatore differintegrale, il differen-ziale a denominatore può essere elevato a potenza senza l’uso di paren-tesi, infatti il simbolo d(x− a) rappresenta un termine infinitesimale enon un prodotto, e non può essere spezzato.

Nella stessa tabella, 1, sono indicati tutti gli altri simboli. Non restache aggiungere che la trasformata L {x(t)} di una funzione x(t), siindica con la stessa lettera della funzione, ma con la maiuscola, cioèX(s).

A C R O N I M I

ODE Ordinary Differential Equation

Equazione Differenziale Ordinaria. Equazione in cui compa-iono esclusivamente integrali e derivate di ordine intero o nullo.In caso contrario l’equazione è detta FODE. Nella definizioneoriginale data in analisi, l’aggettivo ordinaria sta ad indicare ladipendenza della funzione e di tutte le sue derivate, da una solaincognita.

FODE Fractional-Order Differential Equation

Equazione Differenziale di Ordine Frazionario. Un’equa-zione differenziale si dice di ordine non intero se in essa compa-re almeno una derivata o integrale di ordine non intero. In casocontrario l’equazione è detta ODE.

DTFT Discrete Time Fourier Transform

Trasformata di Fourier tempo-discreta. Nota relazione inte-grale per la trasformazione in frequenza di una funzione, svi-luppata dal matematico e fisico francese Jean Baptiste JosephFourier (1768-1830) nel 1822.

ROC Regione di Convergenza

Region Of Convergence. Dominio del piano complesso nelquale la trasformata di Laplace converge ad un valore finito.

FOC Fractional-Order Controller

Controllore di Ordine Frazionario. I controllori fraziona-ri offrono maggiori possibilità di progetto, grazie al grado dilibertà offerto dall’ordine del controllore.

PID Proportional Integral and Derivative

Moltiplicatore Integratore Derivatore. Controllore basilarecostituito da tre blocchi elaborativi, un integratore, un derivatoreed un blocco di guadagno.

PIλDµ Proportional Integral and Derivative

Moltiplicatore Integratore Derivatore (di ordine non inte-ro). Controllore costituito da tre blocchi, analoghi a quelli diun PID classico, ma in cui i blocchi che svolgono le funzione diintegrazione e derivazione hanno ordini qualunque.

TID Tilted Integral and Derivative

Sagomatore Integratore Derivatore. Controllore costituitoda tre blocchi elaborativi, un integratore, un derivatore, ed unblocco per risagomare la caratteristica di Bode. Costituisce unarielaborazione variata del controllore PID.

CRONE Contrôle Robuste d’Ordre Non Entier

Acronimo francese per indicare il Controllo Robusto di Or-dine non Intero, metodologia di controllo sviluppata dal teamguidato dal ricercatore francese Alain Oustaloup [32].

ix

x ACRONIMI

MATLAB © Matlab ©

Matlab © è un programma della MathWorks © Inc, 1994-2009.

LTI Lineari Tempo Invarianti

Classe di sistemi descritti da un’equazione differenziale linearea coefficienti costanti della forma (2.1) o (2.2), che godono delleproprietà di linearità (la funzione e tutte le sue derivate figuranosolo con esponente unitario o nullo) e permanenza (tutti i coeffi-cienti sono costanti). Nei testi in lingua inglese tale proprietà èdefinita come LSI, Linear Shift-invariant.

G-L Grünwald-Letnikov

Tra le definizioni di differintegrali, è ragguardevole l’equazio-ne (1.13) trovata dai matematici Anton Karl Grünwald e Alek-sey Vasilievich Letnikov, che spesso si trova indicata come dif-ferintegrale di G-L.

R-L Riemann-Liouville

Una definizione di differintegrali molto pratica, è l’equazione (1.14)trovata dai matematici Bernhard Riemann e Joseph Liouville,che spesso si trova indicata come differintegrale di R-L.

M-L Mittag-Leffler

Autore della funzione che prende il suo nome, e che può con-siderarsi come una generalizzazione di funzioni esponenzialicomplesse.

IEE Institution of Electrical Engineers

Associazione britannica degli ingegneri elettrici, da non confon-dersi con l’americana IEEE.

SISO Single Input Single Output

Ingresso Singolo Uscita Singola. Sistemi monodimensionali incui si ha un segnale di uscita ed un ingresso.

MIMO Multiple Input Multiple Output

Ingressi Multipli Uscite Multiple. Sistemi a più dimensioniin cui si hanno molteplici segnali di uscita e di ingresso.

FIR Finite Impulse Response

Risposta Impulsiva Finita. Filtri con risposta impulsiva di lun-ghezza finita. Sono semplici da implementare, la stabilità è im-mediata se i coefficienti del filtro sono limitati.

IIR Infinite Impulse Response

Risposta Impulsiva Infinita. Filtri con risposta impulsiva infi-nita, che presentano un meccanismo di feedback per cui l’uscitadipende dall’uscita stessa negli istanti precedenti.

BIBO Bounded Imput Bounded Output

Criterio per la verifica della stabilità dei sistemi. Definisce unsistema stabile, se per ingressi limitati le uscite sono limitate.

SNR Rapporto Segnale Rumore

Signal to Noise Ratio. Rapporto tra la potenza del segnaleutile e la potenza del segnale non desiderato (rumore).

ACRONIMI xi

PSE Power Series Expansion

Espansione in Serie di Potenze. Si intende genericamente unaqualsiasi espansione basata sulle serie di potenze, tipicamen-te l’espansione di MacLaurin, ma in generale l’espansione diTaylor o di Laurent sono altri due esempi di PSE.

CFE Continued Fraction Expansion

Espansione in Frazioni Continue. L’espansione in frazionicontinue è una rappresentazione di un numero reale o di un si-stema frazionario, basato sul successivo annidamento di frazio-ni.

ARMA Auto Regressive Moving Average

Media Mobile Autoregressiva. Filtro per generare segnali ditipo stocastico a partire da un rumore normale opportunamentemodellato dal filtro stesso.

FARMA Fractional Autoregressive Moving Average

Media Mobile Autoregressiva Frazionale. Espansione al ca-so frazionario dei modelli di sistemi di tipo ARMA.

1 B A S I M AT E M AT I C H E

I sistemi trattati in questo testo sono caratterizzati da termini, che neldominio complesso di Laplace, assumono la forma sα, in cui tipicamen-te α non è un numero intero come nei sistemi ordinari, ma appartieneall’insieme dei numeri reali R, o, più raramente, persino all’insiemedei numeri complessi C.

Passando poi nel dominio del tempo, si ottengono equazioni diffe-renziali di ordine non intero (Fractional-Order Differential Equation,FODE): l’ordine di derivazione di tali equazioni è dato infatti dall’espo-nente α della corrispondente equazione nel dominio trasformato.

Le teorie alla base di tali oggetti matematici, nonostante le analogiecon l’analisi infinitesimale ordinaria che vanno ad espandere, sono suf-ficientemente elaborate da richiedere una trattazione specifica. Le basiteoriche della materia sono state poste e consolidate nel corso degliultimi trecento anni; oggi si prestano ad essere applicate in altri ambitidella fisica, della chimica e dell’ingegneria, come vedremo in questotesto!

Cominciamo quindi, partendo proprio dagli strumenti di base delcalcolo farzionario, in questo primo capitolo, principalmente tratto dal[27] salvo ove diversamente specificato. Un diverso approccio dell’ar-gomento, a partire dagli spazi di Hölder, è illustrato da Samko, Kilbase Marichev [36]. Entrambi i testi, attualmente fuori produzione, sonodisponibili nella biblioteca di Matematica della Sapienza.

1.1 FUNZIONE GAMMALa funzione gamma permette di generalizzare la definizione di fatto-riale anche per numeri reali e complessi, non solo interi, e quindi diestendere di conseguenza le formule che fanno uso del fattoriale. An-che i coefficienti binomiali, necessari nelle definizioni generalizzate diderivata, possono essere definiti con la funzione gamma e quindi gene-ralizzati. La definizione risale ad Eulero (Leonhard Euler, 1707–1783), Definizioneed è la seguente[46]:

Γ(x) ≡∫ ∞

0tx−1e−tdt, (1.1)

per x ∈ R+.La funzione gamma completa si ottiene dalle funzioni gamma in-

complete, superiore ed inferiore, rispettivamente:

Γ(a, x) ≡∫ ∞

xta−1e−tdt

γ(a, x) ≡∫ x

0ta−1e−tdt,

dalle quali

Γ(a) ≡ Γ(a, x) + γ(a, x) =∫ ∞

xta−1e−tdt +

∫ x

0ta−1e−tdt,

1

2 BASI MATEMATICHE

-4 -3 -2 -1 0 1 2 3 4 5-10

-8

-6

-4

-2

0

2

4

6

8

10Funzione Gamma

Figura 1: Funzione Gamma Γ(x) per valori reali di x. Per valori di x ∈ N èstato riportato con un cerchietto il valore del fattoriale di x − 1, aconferma visiva della legge di ricorrenza (1.4) Γ(n) = (n− 1)!.

come si verifica facilmente.La formula fattoriale si ottiene con alcuni passaggi, integrando perGeneralizzazione dei

fattoriali parti a partire dalla definizione data (1.1):

Γ(x) =[−tx−1e−t

]∞

0+∫ ∞

0(x− 1)tx−2e−tdt =

= (x− 1)∫ ∞

0tx−2e−tdt = (x− 1)Γ(x− 1),

otteniamo la legge di ricorsione che lega le funzioni Γ di due elemen-ti distanziati di un’unità. Reiterando la legge di ricorsione per unargomento intero n si può scrivere:

Γ(n) = (n− 1)Γ(n− 1) = (n− 1)(n− 2)Γ(n− 2) = · · · == (n− 1)(n− 2)(n− 3) · · · Γ(1) ≡ (n− 1)!, (1.2)

in cui si è osservato che Γ(1) = 1. Dalla definizione (1.1) infatti,

Γ(1) =∫ ∞

0t1−1e−tdt =

∫ ∞

0e−tdt = −

[e−∞ − e0

]= 1.

La formula fattoriale (1.2) ricavata per un intero n, può essere calcolatala funzione Γ(x)equivale al fattoriale! per un reale qualunque x, per cui le relazioni

Γ(n + 1) = n!

Γ(x + 1) = xΓ(x),

(1.3)

(1.4)

ci permettono di estendere la definizione di fattoriale, che si rivelanoestremamente utili in molti casi.

Un’immediata applicazione delle formule fattoriali appena ricavate,è l’estensione della definizione dei coefficienti binomiali, dall’insiemeCoefficienti binomialidei numeri Naturali N, all’insieme dei numeri Reali R o dei numeriComplessi C. Ricordando la definizione dei coefficienti binomiali [45]

Ck ≡(

nk

)≡ n!

(n− k)!k!, (1.5)


e la seguente proprietà [45],

(−1)k(

nk

)=(

nn− k

)=(

k− n− 1n

)(1.6)

che ci permette in sostanza di espandere le definizione dei coefficientibinomiali anche per n ∈ N−, Sostituendo infine il fattoriale con l’au-silio dell’equazione (1.3), la proprietà dei binomiali (1.6) diventa [27]

(−1)k(

nk

)=(

k− n− 1n

)=

Γ(k− n)Γ(−n)Γ(k + 1)

, (1.7)

equazione che si può considerare la definizione di coefficienti bino-miali per n non intero (n ∈ R, n ∈ C), che permette ad esempio diepandere, non solo le potenze intere di un binomio, ma anche quellenon intere (come ad esempio (a + b)π).

Un’altra relazione utile, si ottiene considerando la seguente proprie- Somma di N terminibinomialità della somma di N termini binomiali [45]:

N

∑n=0

(k− n− 1

k

)=(

N − nn

).

Sostituendo ciascun fattoriale con l’equivalente funzione gamma, gra-zie all’equivalenza (1.3), è possibile generalizzare per n non interiottenendo:

N

∑n=0

Γ(n− α)Γ(−α)Γ(n + 1)

=Γ(N − α)

Γ(1− α)Γ(N). (1.8)

Le formule mostrate in questa sezione trovano diretta applicazio-ne in particolar modo nella definizione delle derivate non intere diGrünwald-Letnikov (G-L), come vedremo nelle prossime sezioni ed inparticolare nel § 1.2.2 nella pagina seguente. La funzione gamma è di- Comando gamma(x)

per Matlabsponibile in Matlab con l’istruzione gamma(x), che restituisce il valoreassunto dalla funzione Γ(x)

1.2 DIFFERINTEGRALIIl calcolo differenziale di ordine intero è un pilastro fondamentale dellamatematica oltre che uno strumento essenziale per la fisica, l’ingegne-ria, ed un numero elevatissimo di altre discipline. Sebbene durante ilcorso di studi si affronti relativamente presto, già in età scolastica, ne-gli insegnamenti di ingegneria si considerano esclusivamente derivateed integrali di ordine intero. Tuttavia, possono esistere anche deriva-te ed integrali di ordine non intero, vediamo quindi come trattare conquesti operatori matematici.

1.2.1 Percorso storico

L’idea di derivate di ordine non intero nasce con le derivate stesse. Già l’idea di Liebniz del1965nel 1695 il matematico e filosofo tedesco Gottfried Wilhelm Leibniz

(1646-1716), parla della possibilità di avere derivate di ordine non in-tero, in una lettera indirizzata al suo collega e matematico francese,Guillaume de l’Hôpital (1661-1704).

4 BASI MATEMATICHE

Dall’idea ai primi studi sistematici, si è dovuto attendere quasi undefinizione diLiouville del 1832... secolo e mezzo, prima che il matematico francese Joseph Liouville

(1809-1882) diede l’impulso alla ricerca, trovando, nel 1832, una primadi definizione di derivata frazionaria di una funzione. Liouville consi-derò uno sviluppo in serie di esponenziali di una funzione, e definì laderivata di ordine non intero operando termine a termine, riportando-si quindi al caso intero. L’idea fu ripresa ed ampliata dal matematico...col contributo di

Riemann tedesco Georg Friedrich Bernhard Riemann (1826-1866), che introdus-se una nuova definizione di integrale definito (analoga alla definizionepresentata al § 1.2.3), che era applicabile alle serie di potenze di or-dine non intero. Fu poi Anton Karl Grünwald (Praga 1838-1920) adunificare i loro studi.

Successivamente, lo stesso Grünwald, superò le limitazioni impostedefinizione diGrünwald del 1867 dalla definizione di Liouville, giungendo nel 1867, con la collabora-

zione del matematico russo Aleksey Vasilievich Letnikov (1837-1888),alla definizione più complessa quanto più naturale, ottenuta a partiredal rapporto incrementale, come mostrato nel § 1.2.2, che fu poi ulte-riormente estesa nel 1930 dal matematico americano Emil Leon Post(1897-1954).

Parallelamente sono stati studiati i possibili ambiti di impiego dellateoria che sia andava delineando. Per quanto riguarda l’elettronica, si ègiunti solamente negli ultimi anni a risultati concretamente applicabilie vantaggiosi. In particolare, l’impiego più promettente e sviluppatoè quello dei controllori frazionari, ma questo è argomento di un’altrocapitolo.

1.2.2 Differintegrali di Grünwald-Letnikov

Derivate generalizzate

Partiamo dal caso intero, per trovare una espressione che possa rac-Derivata dallageneralizzazione del

rapportoincrementale

chiudere in se le definizioni di derivate ed integrali. Considerando laben nota definizione di derivata attraverso il rapporto incremetale

d fdx≡ f ′ ≡ lim

δx→0

f (x)− f (x− δx)δx

. (1.9)

Incrementando l’ordine di derivazione si ottiene una formula ricorsi-derivate successiveva che ci permetterà di generalizzare il caso di derivata n-esima Adesempio nel caso n = 2 si ha:

d2 fdx2 ≡ f ′′ ≡ lim

δx→0

f ′(x)− f ′(x− δx)δx

,

in cui sostituendo l’equazione (1.9) della derivata prima si ottiene

d2 fdx2 ≡ lim

δx→0

f (x)− 2 f (x− δx) + f (x− 2δx)δx2

Analogamente per una derivata di ordine 3 si ottiene:

d3 fdx3 ≡ f (3) ≡ lim

δx→0

f (x)− 3 f (x− δx) + 3 f (x− 2δx)− f (x− 3δx)δx3 ,

Si osserva che i coefficienti che compaiono davanti alla funzione se-guono la regola dei coefficienti binomiali con segni alterni, e questo


è ancor più evidente sviluppando manualmente il calcolo delle deri-vate successive; inoltre si osserva che aumentando l’ordine di deriva-zione intervengono punti della funzione a distanze sempre maggioridall’ascissa x considerata. Possiamo dunque scrivere definizione di

derivata n-esima

dn fdxn ≡ f (n) ≡ lim

δx→0

[(δx)−n

n

∑i=0

(−1)i(

ni

)f (x− iδx)

], (1.10)

che è la formula generale di derivazione per ordine intero n, anche senecessita di ulteriori elaborazioni per i nostri scopi. . .

Integrali generalizzati

L’usuale definizione di un integrale, come limite delle somme di Rie- Integrali comesomme di Riemannmann, è data dalla:

d−1 f[d(x− a)]−1

≡ f (−1)(x) ≡

∫ x

af (y)dy ≡ lim

δx→0

[δx

N−1

∑j=0

f (x− iδx)

]

in cui la funzione integrata si ottiene dall’area sottesa dalla funzioneintegranda, calcolata come somma di N rettangoli di base infinitesimaδx ed area δx · f (x). L’aver fissato l’estremo inferiore di integrazio-ne ci permette di eliminare l’indeterminatezza dell’integrale, cioè lafunzione integranda non è definita a meno di una costante.

Analogamente a quanto fatto per le derivate, esplicitiamo il calcolo di ordine superioredegli intergrali di ordine superiore, di secondo grado:

d−2 f[d(x− a)]−2

≡ f (−2)(x) ≡ limδx→0

[(δx)2

N−1

∑i=0

(i + 1) f (x− iδx)

],

poi di terzo grado:

d−3 f[d(x− a)]−3

≡ f (−3)(x) ≡ limδx→0

[(δx)3

N−1

∑i=0

(i + 1)(i + 2)2

f (x− iδx)

],

otteniamo infine l’equazione della definizione di un generico integrale: definizione diintegrale n-esimo

d−n f[d(x− a)]−n

≡ f (−n)(x) ≡ limδx→0

[(δx)n

N−1

∑i=0

(i + n− 1

i

)f (x− iδx)

],

i cui i coefficienti binomiali generalizzano la successione dei coefficien-ti della funzione.

Unificazione delle definizioni

L’equazione così ottenuta, non è uniforme con l’analoga equazione Unificazione delledefinizionidella derivata (1.10), bisogna lavorarci ancora un po’. Nel caso di inte-

grali, dati x, a ed N, si avrà che δx = (x− a)/N. Posso applicare taledefinizione dell’incremento δx anche per le derivate, come se avessicampionato i punti in cui derivare ad intervalli δx. Tra l’altro, restrin-gendo il dominio di derivazione ponendo a come limite inferiore, laderivata esiste comunque.

6 BASI MATEMATICHE

Sostituendo la definizione data dell’incremento δx, otteniamo quindile seguenti formule per derivate di ordine n:

δn fδxn = lim

N→∞

{(x− a

N

)−n N−1

∑i=0

(−1)i(

ni

)f[

x− i(

x− aN

)]}, (1.11)

ed integrali di ordine n:

d−n fd(x− a)−n = lim

N→∞

{(x− a

N

)−n N

∑i=0

(i + n− 1

i

)f[

x− i(

x− aN

)]},

(1.12)

che differiscono solo per l’espressione dei coefficienti binomiali.Ricordando la proprietà dei binomiali (1.6) che qui riscriviamo:

(−1)k(

nk

)=(

nn− k

)=

Γ(k− n)Γ(k− n)Γ(k + 1)

, (1.6)

si osserva come le definizioni (1.11) e (1.12) siano in realtà equivalenti!Si noti inoltre come fino ad ora si sia considerato il caso n intero.

Estensione al caso non intero

La naturale estensione per ordine non intero si ottiene sostituendo leDefinizionedifferintegrali di

Grünwald-Letnikovnotazioni equivalenti dei coefficienti binomiali (1.6), con la più gene-rica espressione indicata in (1.7) che fa uso della funzione Gamma,ottenendo infine:

aDαx f (x) ≡ lim

N→∞

{(x− a

N

)−α 1Γ(−α)

N−1

∑i=0

Γ(i− α)Γ(i + 1)

f[

x− i(

x− aN

)]}

(1.13)

Abbiamo ottenuto un’equazione che unifica le definizioni di derivataed integrale, permette cioè di ottenere entrambi in funzione dell’or-dine di derivazione α ∈ R, che prende il nome di differintegrale diGrünwald-Letnikov (G-L), dai suoi scopritori.

L’equazione (1.13) si può considerare come la più importante tra levantaggi delladefinizione di G-L definizioni di differintegrali, in quanto:

• richiede il minor numero di restrizioni per le funzioni alle qualisi applica;

• evita l’uso esplicito di derivate o integrali della funzione f (x);

• permette di estendere le regole di derivazione al caso non intero;

sebbene questi pregi si paghino al prezzo di una maggiore complessi-tà!

1.2.3 Differintegrali di Riemann-Liouville

La definizione di Riemann-Liouville è una semplice generalizzazio-Integrazionemultipla di Cauchy


ne della formula di integrazione multipla data da Cauchy (Augustin-Louis Cauchy, matematico e ingegnere francese, 1789-1857), che quiricordiamo

d−n fd(x− a)−n =

∫ x

adxn−1

∫ x1

a· · ·

∫ xn−1

af (x0)dx0 =

=1

(n− 1)!

∫ x

a(x− y)n−1 f (y)dy,

che esprime un integrale multiplo come un integrale singolo opportu-namente pesato, tra l’altro con una funzione peso molto semplice!

La formula di integrazione multipla di Cauchy si può facilmente generalizzazione alcaso non interogeneralizzare al caso non intero: il fattoriale si può sostituire con

la funzione gamma (1.3), l’ordine n si generalizza per α qualunque,ottenendo

dα fd(x− a)α

=1

Γ(−α)

∫ x

a(x− y)−α−1 f (y)dy, (1.14)

valida anche per α ∈ C. Tuttavia si deve sempre tenere in mente chela formula è stata ricavata a partire dall’integrazione della funzionef , quindi in tal caso sarà certamente convergente per <{α} < 0. Sidimostra, che per <{α} > 0, la (1.14) è l’inverso della derivata diordine α, la definizione ottenuta dunque è valida per α qualsiasi.

La formula (1.14) fa uso delle differenze (x− y), cioè delle differenze causalità,anti-causalità edefinizione di Weil

all’indietro, per cui è causale. Facendo uso delle differenze in avanti,cioè scrivendo (x + y), si otterrebbe un’analoga definizione anti-causale,che fa uso solo dei campioni in avanti del segnale. Tale definizioneprende il nome di differintegrale di Weil [30].

La definizione (1.14) è la più diffusa, grazie alla sua praticità d’uso vantaggiin molteplici situazioni.

1.2.4 Differintegrali di Courant e Hilbert

La terza definizione di differintegrali considerata è stata proposta nel Riemann-Liuovilleper α = 1/2

1962 dai matematici Courant e Hilbert, ed è la seguente:

d12 f

d(x− a)12≡ 1√

π

ddx

∫ x

a

f (y)dy√x− y

Si può ricavare tale definizione come caso particolare della definizio-ne di Riemann-Liouville, semplicemente ponendo α = 1/2.

1.2.5 Differintegrali di funzioni elementari

Rincuora osservare come alcune derivate fondamentali di ordine interosiano facilmente estensibili anche al caso non intero. Consideriamoalcune semplici funzioni note e calcoliamone la derivata frazionaria.

Funzione costante

Applicando la formula di G-L (1.13) sulla funzione f (x) = 1 si ottiene

aDαx f (x) = lim

N→∞

{(x− a

N

)−α N−1

∑i=0

Γ(i− α)Γ(−α)Γ(i + 1)

},

8 BASI MATEMATICHE

0

0.5

1 −1−0.5

00.5

1

0

0.5

1

1.5

α

Derivate f.ne costante

x

Dα

Figura 2: Differintegrali della funzione costante f (x) = 1, in cui l’estremoinferiore a = 0, il dominio è x ∈ [0, 1], e l’ordine di derivazione α ∈[−1, 1]. In corrispondenza dei valori interi dell’ordine di derivazioneα, sono tracciate le curve di riferimento: per α = 0 ho la funzione nonderivata, per α = 1 ho la derivata di una costante, cioè la funzionenulla, infine per α = −1 ho l’integrale di una costante, cioè una retta.

poi con l’ausilio della (1.8) si semplifica

aDαx f (x) =

(x− a)−α

Γ(1− α),

nella semplice forma di un rapporto!

La funzione costante, per la sua semplicità, si presta molto beneFacilità diinterpretazione del

significato deidifferintegrali . . .

per capire il significato dei differintegrali di ordine non intero. Fis-siamo l’estremo inferiore a = 0, ed osserviamo cosa accade al variaredell’ordine di derivazione α tra [0, 1] grazie al grafico in figura 2

La superficie interpola le linee continue, che corrispondono a deri-. . . che interpolanoderivate ed integrali

ordinarivate o integrali di ordine intero, quindi noti. Per ordine non interosi ha una transizione graduale tra le curve. Per α ∈ [0, 1] le derivatemostrano l’effetto del punto iniziale della funzione, assumendo valorielevati in prossimità di x = 0 (nel grafico sono stati tagliati i valoritroppo alti, delimitati dalla linea curva in verde). Le curve si pieganoin modo sempre più brusco fino ad arrivare ad α = 1, per cui si ha unangolo retto tra l’origine delle x in cui la derivata tende a infinito, ed ipunti successivi, in cui la derivata è nulla.

Potenze

La derivata prima di una potenza è banalmente Dxp = p · xp−1, qual-cosa di molto simile si ottiene nel caso frazionario:

aDαxxp =

Γ(p + 1)Γ(p− α + 1)

· xp−α


Funzione esponenziale

Anche in questo caso se Decx = c · ecx, ottengo qualcosa di simile nelcaso non intero:

aDαxecx = cα · ecx

1.3 PROPRIETÀ DIFFERINTEGRALIPer le derivate di ordine intero sono state sviluppate molteplici pro-prietà, che possono essere generalizzate anche nel caso non intero. For-se la proprietà più importante, la linearità, discende immediatamenteda ciascuna delle definizioni date nel § 1.2 nella pagina 3.

Un’altra formula molto importante nel caso intero esprime la deri-vata di funzioni composte. Nel caso intero si riesce in questo modo,a ricavare l’espressione della derivata di funzioni anche complesse, apartire dalla funzioni elementari note. Nel caso non intero purtroppo,non esiste una formula altrettanto valida, ponendo dei seri limiti allapossibilità di operare su funzioni complesse.

1.3.1 Differintegrali multipli

Cominciamo col vedere come trattare con 2 successive differintegrazio- Studio delcomportamento didue differintegrazionisuccessive

ni di una funzione f (x), partendo dalle identità

dn

dxn

{dN f (x)

dxN

}=

dn+N

dxn+N =dN

dxN

{dn

dxn

},

valide per n, N appartenenti all’insieme dei numeri naturali N. Cer-chiamo il significato di termini analoghi per ordini di derivazione ointegrazione non interi, che ci permettano di operare con differintegra-zioni multiple.

Per funzioni f , che siano N volte differenziabili si può scrivere Partendo dal teoremafondamentale delcalcolo. . .∫ x

af (N)(y)dy = f (N−1)(x)− f (N−1)(a),

dal teorema fondamentale del calcolo.Integrando una seconda volta1, con qualche passaggio:

d−2 f (N)

d(x− a)−2 =d−1

d(x− a)−1

{d−1

d(x− a)−1 f (N)(x)}

=

=d−1

d(x− a)−1

{f (N−1)(x)− f (N−1)(a)

}=

= f (N−2)(x)− f (N−2)(a)−∫ x

af (N−1)(a)dy =

= f (N−2)(x)− f (N−2)(a)− (x− a) f (N−1)(a),

si può generalizzare al caso n ∈N−, . . . otteniamo la 1a

soluzione, n integraliseguiti da Nderivate;

1 Si ricorda che la notazione delle derivate (e quindi degli integrali) espressa come rap-porto dei differenziali della funzione è dovuta a Leibniz. Qui viene convenientementeusata per esprimere in modo compatto integrali multipli.

10 BASI MATEMATICHE

d−n f (N)

d(x− a)−n = f (N−n)(x)−n−1

∑k=0

(x− a)k

k!f (N+k−n)(a). (1.15)

Si osservi che per N = n la sommatoria non è nulla, a meno che nelpunto x = a non sia nulla la funzione f con le prime n − 1 deriva-te, cioè derivata ed integrale non sono operatori esattamente inversi.Come noto infatti, rende nulli i termini costanti, l’integrale non puòrecuperare il valore perduto dalle operazioni precedenti.

Ripartiamo dalla (1.15) fissando N = 0,derivando la formulaottenuta, ponendo a

0 il numero diderivazioni. . .

d−n fd(x− a)−n = f (−n)(x)−

n−1

∑k=0

(x− a)k

k!f (k−n)(a), (1.16)

e deriviamo l’equazione semplificata così ottenuta. Derivando la primavolta si ha

d f (−n)

d(x− a)= f (1−n)(x)−

n−1

∑k=1

(x− a)k−1

(k− 1)!f (k−n)(a),

infatti il primo termine è integrato una volta in meno, cioè 1− n, laderivata dei termini della sommatoria è la derivata rispetto ad x di untermine (x− a)k.

Continuando a derivare, fino ad N volte, si ottiene. . . otteniamo la 2a

soluzione, N derivateseguite da n integrali

;dN f (−n)

d(x− a)N = f (N−n)(x)−n−1

∑k=N

(x− a)k−N

(k− N)!f (k−n)(a), (1.17)

in cui la sommatoria scompare e resta

dN f (−n)

d(x− a)N = f (N−n)(x), (1.18)

se l’ordine di derivazione N è maggiore dell’ordine di integrazione,N ≥ n.

Ripartendo infine dalla (1.16), e sostituendo −n = N − n otteniamoderivando infineN − n volte . . . la formula composta

dN−n fd(x− a)N−n = f (N−n)(x)−

n−N−1

∑k=0

(x− a)k

k!f (k+N−n)(a), (1.19)

che finalmente chiude il cerchio, ponendo l’ultimo tassello che comple-ta le teorie e risponde alle nostre domande.

Infatti, se nella (1.19) sostituiamo nella sommatoria k con k− N, la. . . si ottengono 3formule equivalenti,

sebbene l’ordine delleoperazioni sia

differente

(1.19) è identica alla (1.17). Non solo, la (1.19) risulta anche essereidentica alla (1.15), per cui, riassumendo, si ha2:

dN f (−n)

d(x− a)N =dN−n f

d(x− a)N−n =d−n f (−N)

d(x− a)−n =

= f (N−n)(x)−n−1

∑k=N−n

(x− a)k

k!f (k+N−n)(a),

relazioni valide nel caso intero.Si osservi che se il limite superiore della sommatoria viene aumen-

tato, compaiono degli integrali calcolati nel punto x = a, tali inte-grali non comportano alcuna variazione dei risultati ottenuti. Questo


ci permette di generalizzare al caso non intero le formule trovate, in Caso non interoparticolare il caso più utile si ha per n = α, con α ∈ R. Si ottiene:

dN f (α)

d(x− a)N =dN+α f

d(x− a)N+α=

dα f (N)

d(x− a)α=

= f (N+α)(x)−N−1

∑k=0

(x− a)k−α−N

Γ(k− α− N + 1)f (k)(a). (1.20)

1.3.2 Differintegrali del prodotto di funzioni

Un noto risultato, del calcolo differenziale di ordine intero, è dato dalla Formula diLeibniz. . .formula di Leibniz per la derivata del prodotto di due funzioni, secondo

la quale

dn

dxn { f · g} =n

∑i=0

(ni

)dn−i fdxn−i ·

digdxi . (1.21)

Anche in questo caso si osserva che il risultato non cambia se la . . . e generalizzazioneal caso non interosomma viene estesa fino ad infinito, dal momento che i coefficienti

binomiali (ni ) sono nulli per i > n. Per tanto, l’estensione naturale

della formula di Leibniz si dimostra essere

dα

d(x− a)α { f · g} =n

∑i=0

(α

i

)dα−i f

d(x− a)α−i ·dig

d(x− a)i , (1.22)

in cui n è il primo intero più grande dell’ordine di derivazione α, cioèn− 1 < α < n.

1.3.3 Differintegrali di funzioni composte

La nota regola per le funzioni composte di ordine intero La formula per lederivate composte . . .

ddx{

f(

g(x))}

=d

dg(x)f(

g(x)) d

dxf (x), (1.23)

è priva di una semplice controparte nel calcolo integrale.Dal momento che nei differintegrali l’ordine può essere tanto posi- . . . non ha semplici

equivalenti per gliintegrali,

tivo quanto negativo, si intuisce che, una generalizzazione sufficiente-mente semplice da poter essere facilmente sfruttabile in pratica, nonesiste.

La generalizzazione per derivate di ordine n si ottiene nel modo per cui lageneralizzazione alcaso n qualunque. . .

usuale, derivando la (1.23) più volte, fino a trovare una formula gene-rale con l’ausilio della formula di Leibniz (1.21). Derivando la primavolta

ddx

{d f(

g(x))

dx

}=

ddx

{d f(

g(x))

dg(x)d

dxg(x)

}=,

2 NB: La formula è presa dal [27], nel quale manca il termine f N−n(x) qui aggiunto percoerenza.

12 BASI MATEMATICHE

poi semplificando la scrittura ponendo f(

g(x))

= f (u) ,

=d

dx

{d f (u)

dudg(x)

dx

}=

ddx{

f ′(u) g′(x)}

=

= g′(x)d f ′(u)

dx+ f ′(u)

dg′(x)dx

=

= g′(x) · d f ′(u)du

dg(x)dx

+ f ′(u)g′′(x) =

= f ′(u)[g′(x)

]2 + f ′(u)g′′(x),

in cui il termine d f ′(u)/dx è stato a sua volta esplicitato con la (1.23).Ripetendo il procedimento per ordini successivi, si ottiene la formu-. . . con la formula di

Faà di Bruno. . . la di Faà di Bruno per la generalizzazione al caso n-esimo

dn f(

g(x))

dxn = n!n

∑m=1

f (m) ∑P

n

∏k=1

1Pk!

(g(k)

k!

)Pk

,

in cui la somma su P si estende su tutte le possibili combinazioni di ninteri non negativi P1, P2, . . . , Pn, che soddisfano le condizioni{

∑nk=1 kPk = n,

∑nk=1 Pk = m.

Per complicare ancora un po’ la formula di Faà di Bruno, generaliz-. . . per ordinequalunque. . . ziamo al caso non intero [27, pagina 80]

dα f(

g(x))

d(x− a)α=

(x− a)−α

Γ(1− α)f(

g(x))+

+∞

∑i=1

(α

i

)(x− a)i−α

Γ(i− α + 1)i!

i

∑m=1

f (m) ∑P

i

∏k=1

1Pk!

(g(k)

k!

)Pk

,

in cui compare un addendo iniziale ed una sommatoria aggiuntiva diinfiniti termini.

La difficoltà pratica dell’utilizzo della formula trovata di fatto ne. . . non permetteapplicazioni pratiche

per l’eccessivacomplessità!

impedisce qualsiasi reale applicazione. Di fatto ciò costituisce un osta-colo enorme alla differintegrazione di funzioni complesse che non sia-no direttamente ricavabili, ad esempio attraverso la definizione di G-L.Infatti, uno degli strumenti più potenti per la derivazione nel casointero, è proprio dato dalla proprietà (1.23), l’assenza di uno strumen-to ad essa equivalente è una grave carenza della nuova teoria delladifferintegrazione generica!

1.4 TRASFORMATA DI LAPLACELa trasformata di Laplace3

L {x(t)} ≡ X(t) =∫ ∞

0x(t)estdt

è uno strumento indispensabile per lo studio dei sistemi elettronici. LoI sistemi descritti daequazioni

differintegrali simodellano con

sistemi in s conesponenti non interi!

3 Nel 1785 Pierre-Simon Laplace (1749–1827), matematico fisico e astronomo francese,riprendendo gli studi iniziati nel 1744 da Eulero sulle trasformazioni integrali, diede ladefinizione di trasformata che conosciamo.

1.4 TRASFORMATA DI LAPLACE 13

scopo di questa sezione è capirne il funzionamento per ordini non inte-ri. Fortunatamente si possono ottenere risultati del tutto simili al casointero: derivate frazionarie si trasformano in s in potenze frazionariedi grado α, come mostrato nella (1.26). Come diretta conseguenza iteoremi del valore iniziale e finale non cambiano. Vediamo in questaparte come si modificano le proprietà fondamentali delle trasformatedi Laplace.

1.4.1 Convoluzione e prodotto

Anche nel caso frazionario, si ha la corrispondenza tra la convoluzione La convoluzione neltempo si trasformanel prodotto, ancheper sistemi noninteri!

nel tempo ed il prodotto in frequenza [24]. La dimostrazione noncambia per esponenti non interi, pertanto non sarà qui ripetuta, siriporta solo il risultato raggiunto:∫

x(τ)y(t− τ)dτ ≡ x(t)⊗ y(t)⇔ L {x(t)}L {y(t)} ≡ X(s)Y(s).

(1.24)

1.4.2 Proprietà della derivata

Cerchiamo la generalizzazione della proprietà della trasformata delladerivata, cercando di capire quindi cosa si ottiene dalla trasformata diLaplace di un differintegrale.

Indichiamo con L la trasformata di Laplace, x(t) una funzione neldominio del tempo e X(s) la sua corrispondente trasformata.

La nota relazione di ordine intero Analogia tra letrasformate di ordineintero e non

L

{dnx(t)

dtn

}= snX(s)−

n−1

∑k=0

sn−1−kDkx(t)∣∣t=0 , (1.25)

in cui la sommatoria scompare nel caso di integrali, si estende al casonon intero

L

{dαx(t)

dtα

}= sαX(s)−

n−1

∑k=0

skDα−1−kx(t)∣∣t=0 , (1.26)

in cui n è il primo intero maggiore di α, per cui n− 1 ≤ α ≤ n. I terminidella sommatoria rappresentano le condizioni iniziali del sistema, sesi considera un sistema in condizioni iniziali di riposo, tutti i terminisaranno nulli. Vediamo come ricavare la (1.26).

Partiamo considerando indici α negativi, per cui possiamo scrivere Dimostrazione perα < 0: Riemann-Liouville. . .

per la definizione di Riemann-Liouville (1.14):

dα fd(x− a)α

=1

Γ(−α)

∫ x

a(x− y)−α−1 f (y)dy, per q < 0. (1.14)

L’integrale può essere visto come la convoluzione (1.24) tra la funzione . . . con la proprietà diconvoluzione;f (x) e la funzione x−α−1. In tal caso, la trasformata di Laplace della

(1.14), sarà pari al prodotto delle trasformate delle funzioni,

L

{∫ x

a(x− y)−α−1 f (y)dy

}=

= L{

x−α−1 ⊗ f (x)}

= L{

x−α−1}

L { f (x)} (1.27)

14 BASI MATEMATICHE

per cui ricordando che

L{

x−α−1}

=(−α− 1)!

s−α,

si ha

L

{dα f

d(x− a)α

}=

1Γ−α

L{

x−α−1}

L { f (x)} = sαL { f (x)} , (1.28)

in cui si è fatto uso della proprietà (1.4), e la derivata è stata espres-sa con la (1.14) con la trasformata dell’integrale modificato secondola (1.27). Abbiamo dunque ottenuto per la trasformata di un integralefrazionario, una formula identica alla (1.25) valida nel caso intero.

Passiamo adesso al caso più generico per α > 0, in tal caso alla defini-Dimostrazione perα > 0:

Riemann-Liouville digrado n + (α− n). . .

zione di Riemann-Liouville occorre specificare la condizione della leg-ge di composizione (1.20), per cui è possibile sostituire α→ n +(α− n),spezzando poi banalmente le derivate

dα f (x)dxα

=dn

dxn

{dα−n f (x)

dxα−n

}=

dn

dxn {ϕ(x)} , (1.29)

possiamo considerare la differintegrazione di ordine α− n di f comeuna funzione autonoma, in questo modo non resta che la derivata diordine intero.

In tal caso, la trasformata di Laplace della (1.29) si calcola nel modo. . . in parteequivalente al caso

α < 0. . .usuale con la (1.25), per cui

L

{dα f (x)

dxα

}= L

{dn

dxn {ϕ(x)}}

=

= snL {ϕ(x)} −n−1

∑k=0

sk dn−1−k

dxn−1−k

{dα−n f (x)

dxα−n

}. (1.30)

Il primo termine è la trasformata di Laplace di un termine derivatoα− n volte, poiché stiamo considerando 0 < n− 1 < α < n, il termineα− n sarà negativo, per cui

snL

{dα−n f (x)

dxα−n

}= sn {sα−nL { f (x)}

}= sαF(s), (1.31)

come risulta applicando la (1.28).Per ciascun termine della sommatoria invece, sfruttiamo la formula-. . . in parte con la

legge di composizione. . .

zione esplicita delle derivate successive (1.20), in cui si identifica:

N → n− 1− k,

α→ α− n,

}⇒ N + α → α− 1− k,

N−1

∑h=0→

n−k−2

∑h=0

,

in particolare, per l’indice h della sommatoria presente nella (1.20), conk = n− 1 può assumere al massimo il valore n− k− 2 = −1, per cuila sommatoria scompare e della (1.20) non resta che il primo addendo:

dn−1−k

dxn−1−k

{dα−n f (x)

dxα−n

}=

dα−1−k f (x)dxα−1−k . (1.32)


Rimettendo infine insieme tutti i pezzi (1.31), (1.32) nella (1.30) si . . . fino a dimostrarela tesi!ottiene

L

{dα f (x)

dtα

}= sαF(s)−

n−1

∑k=0

skDα−1−k f (x)∣∣t=0 (1.26)

come ci eravamo proposti di dimostrare.

1.4.3 Teorema del valore iniziale e finale

La proprietà della derivata (1.26), ci permette di ricavare il teorema delvalore iniziale, e del valore finale [28], con una dimostrazione identicaalla stessa introdotta nel caso intero.

Teorema del valore iniziale

Per il teorema del valore iniziale,

limt→0 aD

α−1t x(t) = lim

s→∞sαX(s), <(s) > 0.

Teorema del valore finale

Analogamente, anche il teorema del valore finale,

limt→∞ aD

α−1t x(t) = lim

s→0sαX(s), <(s) > 0,

ha una formulazione del tutto identica per sistemi di ordine α reale ointero.

1.5 SUPERFICI DI RIEMANNAbbiamo visto nel § 1.4 nella pagina 12, che i sistemi di ordine nonintero si modellano con funzioni nel dominio di Laplace con esponentenon intero. In tal caso, le potenze diventano funzioni a più valori sulpiano complesso: per un solo termine sα si hanno più punti del pianocomplesso che soddisfano l’equazione.

Le superfici di Riemann [39] permettono di trattare funzioni a più Superfici di Riemann:funzioni a più valoricome se fossero asingolo valore

valori, come se fossero funzioni a singolo valore, con evidenti vantaggiin termini di semplificazioni dei formalismi matematici.

Si ricorda che una funzione si dice a più valori (o polidroma o polivo-ca) se per un punto nel dominio corrispondono più valori distinti nelcodominio [39]. Tra le funzioni a più valori si hanno il logaritmo e laradice.

1.5.1 Radice

Nel campo complesso, dato un un qualsiasi elemento z ∈ C è possi- La fase della radicen-esima si ottiene conn formule distinte

bile trovare esattamente n sue radici distinte, con fase equi-spaziatae modulo identico. Con la notazione polare z = rejϕ, le n radiciwk = n

√z = rejϑk , con k = 0, 1, . . . , n− 1, sono date da:{r = n√

rϑk = ϕ

n + 2kπn , k = 0, 1, . . . , n− 1,

(1.33)

16 BASI MATEMATICHE

come noto dai corsi di analisi.Chiariamo meglio il concetto, considerando, senza perdita di gene-

ralità, le 2 radici per n = 2. Queste, come si ottiene dalle (1.33), sonocaratterizzate da modulo r =

√r, e fase

ϑ =

{ϕ/2 per k = 0ϕ/2 + π per k = 1

(1.34)

Considerando come dominio l’intero piano complesso, i punti z ∈ C

possono essere considerati con fase −π < arg(z) ≤ π, ma il valoreprincipale della radice w0 =

√z =

√rejϑ1 , ϑ1 = ϕ/2, copre solo il

semipiano a parte reale positiva C+. Affinché anche il codominio co-pra l’intero piano complesso, occorre considerare anche i punti confase ϑ2 = ϕ/2 + π, ma in tal caso le soluzioni sono date da 2 differentiformule.

L’insieme di soluzioni di ciascuna delle formule risolventi, è dettorami e linee di taglioramo della funzione, per cui la funzione radice quadrata avrà 2 rami.La definizione della fase in un intervallo limitato, suppone che sianopresenti delle linee di confine coincidenti con gli estremi dell’intervallodelle fasi considerato. Tali linee sono dette linee di taglio, nel casoconsiderato si è scelto, come linea di taglio, la semiretta dei numerireali negativi.

Un modo per considerare la radice quadrata come funzione ad unsolo valore, è di duplicare il dominio e di rimuovere le discontinuità.Per −2π < arg(z) ≤ 2π ho considerato come dominio un sovrapposi-zione di due piani complessi, in tal caso il codominio della funzioneconsiderata, la radice quadrata, coincide con l’intero piano complessoC, pur considerando il solo valore principale, quindi una sola formula.

L’intuizione di Riemann è stata proprio quella di definire un do-Riemann: con ndomini

opportunamentecollegati, calcolo

tutte le radici conuna sola formula

minio in modo opportuno in modo da rendere univoca e continuala soluzione di una funzione. Consideriamo due piani complessi C

sovrapposti (o fogli), e supponiamo di tagliare il primo, in modo dapoter passare, dopo una rotazione completa in senso antiorario attornoall’origine (1° ramo), al piano complesso successivo ad esso sovrappo-sto (2° ramo), e così di nuovo, dopo una ulteriore rotazione di 2π dipoter tornare al punto di partenza. Tale superficie è detta superficie diRiemann della radice quadrata.

Se si considera come dominio la superficie di Riemann, il valore prin-cipale della radice quadrata ha come codominio il piano complesso,che copre con continuità e senza ambiguità!

L’ambiguità residua del metodo di Riemann è data dall’arbitrarietàCiascuna funzioneammette infinite

superfici di Riemann!della scelta delle linee di taglio. Infatti, se consideriamo sempre il casodella radice, esistono infinite linee di taglio che permettono di costruireinfinite superfici di Riemann, date dalle semirette uscenti dall’origine.La superficie principale di Riemann, è data da una scelta particolare dellelinee di taglio, che nel caso della radice corrisponde al semi-asse realenegativo [28].

1.5.2 Logaritmo

Nel campo complesso, anche il logaritmo è una funzione a più valori,Il logaritmo in C èuna funzione ad

infiniti valoriper cui è possibile sfruttare le superfici di Riemann. Si ricorda che datoun numero complesso z, il suo logaritmo è definito come segue:

log(z) = ln |z|+ j arg(z), (1.35)


−1−0.5

00.5

1 −1−0.5

00.5

1−1

−0.5

0

0.5

1

ℑ(z)

=

(z)

ℜ(z)

ℜ(w

)

Superficie di Riemann per f

−1

0

0.5

1

ℑ(w)

−0.5

Figura 3: Superficie di Riemann per la radice w di un numero complesso zottenuta in Matlab. Il dominio è dato dalle z tali che−2π < arg(z) ≤2π. Sull’asse z la parte reale di w, mentre in colore la sua parteimmaginaria.

in cui l’argomento è definito a meno di multipli di 2π, pertanto è unafunzione ad infiniti valori.

La superficie di Riemann in questo caso, si costruisce con infiniteripetizioni del procedimento adottato nel caso della radice quadrata.Consideriamo infiniti fogli sovrapposti, tagliamo ciascun foglio in cor-rispondenza dell’asse dei numeri reali positivi, poi colleghiamo in mo-do continuo ciascun foglio, cioè per fase crescente si passa ai fogli so-vrapposti collegando le rispettive linee di taglio. Si può intuire che inquesto modo otteniamo una spirale che si avvolge come la filettaturadi una vite.

1.5.3 Rappresentazione grafica

È possibile rappresentare le superfici di Riemann graficamente. La rap-presentazione ideale dovrebbe mostrare la corrispondenza tra i puntidella superficie di Riemann ed i punti dell’immagine della funzionecalcolata. È evidente che ciò richiederebbe almeno 4 dimensioni peruna corrispondenza da un piano complesso (a 2 dimensioni) ad unaltro piano complesso.

Potendo visualizzare solo grafici al massimo con 3 dimensioni, si Superfici di Riemannrappresentate comefogli sovrapposti,distanziati delmodulo dellafunzione e con scaladi colori della fase

può indicare con una scala di colori una delle dimensioni coinvolte.Una tipica rappresentazione che rende bene l’idea della superficie diRiemann, riferendosi al caso della radice, è di tracciare un grafico, incui si fa corrispondere a ciascun punto z del piano complesso, la partereale della sua radice <{w} = <{

√z}, mentre la parte immaginaria è

resa con una scala di colore.Ne è un esempio il grafico 3, che rende visivamente il doppio piano

complesso sovrapposto, con continuità della fase tra −2π e 2π. Il co-dice che la realizza è mostrato nel listato, preso dal comando Matlab

cplxroot:

1 r = (0:m)’/m; % Modulo di z

18 BASI MATEMATICHE

2 theta = pi*(-n*m:n*m)/m;

3 z = r * exp(i*theta);

4 s = r.^(1/n) * exp(i*theta/n);

5

6 surf(real(z),imag(z),real(s),imag(s));Listato 1.1: Comandi Matlab per generare la superficie di fig. 3

in cui m è il numero di punti radiali considerati (nell’esempio m =20), n è l’ordine della radice (per la radice quadrata n = 2). È possibilegenerare con gli stessi comandi, le superfici di Riemann delle radici diordine superiore, agendo sulla variabile n.

È possibile inoltre tracciare con la stessa logica il grafico del logarit-mo, che è dato da un piano a forma di spirale infinita.

1.6 FUNZIONE E DI MITTAG-LEFFLERUna funzione, che scopriremo presto essere di estrema importanza,è la funzione definita dal matematico svedese Gösta Mittag-Leffler(1846-1927) che prende il suo nome. Può essere assimilata ad unageneralizzazione della funzione esponenziale [30], o più in generaledi tutte le funzioni esprimibili come combinazioni lineari di funzio-ni esponenziali complesse, per questo è usata nell’ambito dei sistemifrazionari, come vedremo in seguito.

La definizione a due parametri della funzione di Mittag-Leffler (M-L),Definizioneindicata con il simbolo Eα,β, è data dalla seguente serie [47]:

Eα,β(z) ≡∞

∑k=0

zk

Γ(αk + β),

valida per α > 0, β > 0.Spesso la funzione di M-L è indicata con un solo parametro, in talanalogia della

funzione di M-L conl’esponenziale e, da

cui deriva ilsimbolo E

caso si intende β = 1. Se anche α = 1, ricordando che k è interoe ricordando la formula (1.3) di ricorsione della funzione Gamma, siottiene:

E(z) ≡∞

∑k=0

zk

Γ(k + 1)=

∞

∑k=0

zk

k!≡ ez, (1.36)

espressione che coincide con la definizione della funzione esponenzia-le, come mostrato nell’ultimo passaggio.

Recentemente, è stata proposta dal prof. Podlubny [34] un’imple-mentazione efficiente in Matlab della funzione di M-L, disponibile sulsito Matlab Central

4. Le curve in figura 4 nella pagina successivamostrano l’andamento della funzione, ottenuto con lo script citato.

1.6.1 Trasformata di Laplace di E

Operando sui sistemi LTI (lineari tempo-invarianti) nel dominio di La-Cerchiamo lacorrispondenza tratempo e frequenza,

per sistemi di ordinenon intero

place, le funzioni che si ottengono sono rapporti di polinomi in s. Ènoto dall’analisi di sistemi interi che sussiste la corrispondenza

L{

tke−at}

=(k− 1)!(s + a)k ,

4 La funzione mlf è disponibile all’indirizzo http://www.mathworks.com/matlabcentral/fileexchange/8738.

http://www.mathworks.com/matlabcentral/fileexchange/8738


1.6 FUNZIONE E DI MITTAG-LEFFLER 19

−30 −20 −10 0 10−1.5

−1

−0.5

0

0.5

1

1.5

2

2.5Funzione di Mittag-Leffler

a=1, b=1a=2, b=1a=1, b=4a=3, b=1

Figura 4: Grafico della funzione di M-L tracciato con lo script mlf, per alcunivalori dei parametri α e β.

che ci permette di passare rapidamente dal dominio del tempo al do-minio trasformato e viceversa, senza ricorre alla definizione esplicitadi trasformata e di antitrasformata.

Ci proponiamo di ricavare un’analoga relazione, che possa essereapplicata anche per sistemi di ordine non intero. Una spiegazionesintetica per β = 1 è disponibile su [47], seguiamo invece l’ampia ecompleta dimostrazione data da [34].

La trasformata di Laplace di un termine esponenziale semplice e−t, Procedimento:trasformiamo e−t,con lo sviluppo inseriedell’esponenziale. . .

come noto, coincide con un termine razionale in s. Partiamo proprioconsiderando la seguente funzione integrale,∫ ∞

0e−te±ztdt =

∞

∑k=0

±zk

k!

∫ ∞

0e−ttkdt =, (1.37)

che ci proponiamo di risolvere con la definizione della funzione espo-nenziale (1.36), in modo differente da quanto fatto nei corsi base diingegneria. Si osservi come si sia utilizzato il simbolo z al posto delconsueto s, per indicare la variabile di Laplace.

L’integrale (1.37) può essere risolto per separazione di variabili o perassociazione ad una trasformata di Laplace nota: considerando l’argo-

mento e−sttkdt, posso sfruttare la nota formula L{

tk}

= k!/sk+1, checalcolata per s = 0 e sostituito nell’espressione (1.37) si ottiene

∞

∑k=0

±zk

k!k! =

∞

∑k=0

(±z)k =1

1∓ z, (1.38)

sfruttando i risultati noti per la somma delle serie geometriche. Riepi- . . . otteniamo la notaequazione (1.39) o(1.40). . .

logando si ha l’uguaglianza∫ ∞

0e−te±ztdt =

11∓ z

. (1.39)

La stessa uguaglianza con la funzione d M-L sarà∫ ∞

0e−tEα,β(±zt)dt =

11∓ z

, (1.40)

20 BASI MATEMATICHE

in cui chiaramente α, β = 1.Deriviamo ora k volte la (1.39), ottenendo. . . che derivata k

volte. . . ∫ ∞

0e−ttke±ztdt =

k!(1∓ z)k+1 , (1.41)

poi applichiamo gli stessi risultati all’equazione (1.39) che fa uso dellafunzione di M-L,∫ ∞

0e−ttβ−1Eα,β(±ztα)dt =

k!(1∓ z)k+1 ,

notando sempre che k rappresenta l’ordine di derivazione.Confrontiamo infine entrambe le equazioni, con la trasformata note-

vole:

L{

tke±at}

=∫ ∞

0e−sttke±ztdt =

k!(s∓ z)k+1 ,

ottenendo infine. . . fornisce latrasformata di

Laplace completadella funzione di

M-L!

L{

tαk+β−1E(k)α,β(±tαa)

}=∫ ∞

0tαk+β−1E(k)

α,β(±tαa)e−stdt =k!sα−β

(sα ∓ a)k+1

.

Tale equazione ci permette di antitrasformare un generico termine del-la forma del rapporto di polinomi in s di esponente non intero, edesprimerne la sua antitrasformata in termini di derivate della funzio-ne di M-L! In particolare per antitrasformare un termine razionale ins, ho:

L −1{

k!sα−β

(sα ∓ a)k+1

}= tαk+β−1E(k)

α,β(±tαa), (1.42)

equazione semplificabile come spesso accade quando il termine in sal numeratore non è presente, in tal caso α = β, e l’esponente deldenominatore è pari ad 1, cioè k = 0 ottenendo la seguente formula diSpesso è sufficiente la

versione semplificata antitrasformazione

L −1{

1sα ∓ a

}= tα−1Eα,α(±tαa), (1.43)

di estrema utilità!

2 S I S T E M I A N A LO G I C I

Dopo l’introduzione degli strumenti matematici di base vista nel pri-mo capitolo, cominciamo lo studio di sistemi non interi nel dominiocontinuo del tempo. Sebbene si abbia a che fare con concetti e strumen-ti di calcolo simili agli stessi sviluppati per i sistemi interi, la maggioregeneralità conquistata si paga a prezzo di una maggiore complessitàanalitica: ad esempio occorre trattare con funzioni a valori multipli, ocon la funzione di Mittag-Leffler al posto dell’esponenziale.

2.1 DEFINIZIONI GENERALIUn sistema analogico lineare tempo invariante può essere modellato,nel dominio temporale, con un’equazione lineare alle differenze finite,della forma

N

∑n=0

anDny(t) =M

∑m=0

bmDmx(t), (2.1)

in cui i coefficienti an e bn sono costanti fissate e indipendenti dal tempo.I circuiti rappresentabili dall’equazione (2.1) sono detti Lineari TempoInvarianti (LTI).

Considerando lo stesso sistema per derivate non intere, siano αn gli Tra lerappresentazioni deisistemi interi e non,si hanno fortianalogie!

N ordini di derivazione, la scrittura

N

∑n=0

anDαn y(t) =M

∑m=0

bmDαm x(t), (2.2)

generalizza la (2.1), per αn ∈ R.Applicando la proprietà (1.26) introdotta nel § 1.4 nella pagina 12,

per ciascuno dei termini della sommatoria, si ottiene

N

∑n=0

ansαn Y(s) =M

∑m=0

bmsαm X(s), (2.3)

valida per circuiti con condizioni iniziali di riposo, cioè per x(t)∣∣t=0 =

0. L’equazione ottenuta viene tipicamente riscritta sotto forma di rap-porto di polinomi, vale a dire come funzione di trasferimento: anche per la funzione

di trasferimento

H(s) = ∑Mm=0 bmsαm

∑Nn=0 ansαn

. (2.4)

La stessa espressione permette di ottenere la risposta in frequenza conla sostituzione ejω → s. Per tracciare i diagrammi ideati dall’ingegne-re americano Hendrik Wade Bode (1905-1982), occorre calcolare comedi consueto, il modulo e la fase della risposta in frequenza H(jω). Ilmodulo è costante per tutte le radici di jω, mentre per la fase tipi-camente si considera il solo valore principale, o equivalentemente la

21

22 SISTEMI ANALOGICI

prima intersezione della superficie di Riemann con l’asse dei numeriimmaginari positivi.

Osserviamo come nel caso non intero la pendenza del modulo suldiagramma di Bode, possa assumere pendenze qualsiasi. Infatti, se nelI sistemi non interi

possono assumerequalsiasi pendenzasul diagramma di

Bode!

caso intero a distanza di una decade il modulo varia di un multiplo di20 dB, infatti per la proprietà della somma dei logaritmi

20 log10((10ω)n) = 20n log10(ω)+ 20n log10(10) = 20n log10(ω)+ 20n,

si osserva che a distanza di una decade mi sposto di 20n dal valo-re precedente. Se al posto di n ho α ∈ R, posso avere un valorequalunque!

Per lavorare con i sistemi non interi, è necessario disporre degli stessistrumenti tipicamente utilizzati con i sistemi di ordine intero, bisognaquindi capire come svolgere le seguenti operazioni:

• Calcolare la funzione di trasferimento. In realtà basta trasforma-re l’equazione differenziale (2.1) come fatto nell’introduzione perricavare la (2.4);

• Antitrasformare la funzione di trasferimento o l’uscita del siste-ma. Trovata l’uscita risolvendo linearmente l’equazione nel do-minio trasformato, bisogna tornare nel tempo, per avere la ri-sposta impulsiva o la y(t). È possibile farlo tramite scomposi-zione in fratti semplici ed antitrasformate elementari (una voltacalcolate!), così come si fa per sistemi di ordine intero;

• Tracciare i diagrammi di Bode;

• Trovare la risposta forzata per condizioni iniziali non nulle;

• Studiare la stabilità dei sistemi.

nelle prossime pagine capiremo come si eseguono le operazioni de-scritte, anche per sistemi di ordine non intero.

2.2 SISTEMI ELEMENTARIAnche per i sistemi di ordine non interi, il modo più pratico di pro-cedere consiste nella scomposizione in sistemi elementari, con un pro-cedimento del tutto simile alla scomposizione in frazioni parziali tipi-camente usata per i sistemi di ordine intero. Si prende un sistema, losi scompone in una somma di sistemi elementari, poi si antitrasformaciascuno di questi per ottenere la risposta nel dominio del tempo, orisposta impulsiva del sistema.

Per prima cosa vediamo quindi le caratteristiche dei sistemi ele-mentari non interi, poi cerchiamo di capire come invertire un sistemaqualunque, dal dominio complesso di Laplace, al tempo t.

2.2.1 Termine monomio

Il sistema più frazionario più semplice a cui possiamo pensare, ha unafunzione di trasferimento caratteristica

H(s) = sα, α 6= 0, (2.5)


ricordando però che l’equazione di trasferimento ha senso e risultadefinita se all’equazione del sistema viene associata una regione diconvergenza!

Una risposta del genere è riscontrabile, ad esempio, per un conden- Frattanza, cioècondensatore oinduttore frazionario

satore o per un induttore di ordine non intero; un dispositivo di questotipo viene detto frattanza, parola che nasce dall’incrocio di induttanzae frazionario. Poi, siccome per grado α qualunque, la distinzione ma-tematica dell’espressione della capacità ed dell’induttanza è solo unadiscriminazione del segno dell’esponente α, il solo termine frattanzava bene per indicare entrambi i dispositivi. La realizzazione fisica diun simile dispositivo sarà vista nel capitolo 4.

Il sistema ha un polo, oppure uno zero nell’origine, in base al segnodell’esponente α.

Come visto nel § 1.5 nella pagina 15, esistono più valori di s che perun H(s) fissato soddisfano la relazione sα = H(s). Alla (2.5) posso-no essere associate infinite superfici di Riemann, per cui l’equazionerappresenta infiniti sistemi, ma un sistema reale può giacere solo sullasuperficie principale di Riemann.

Risposta impulsiva

Per riportare la (2.5) nel dominio del tempo, partiamo dalla definizione Partendo dalladefinizione diderivatafrazionaria. . .

integrale di Riemann-Liouville (1.14) data nel § 1.2.3 nella pagina 6,valida per α < 0

aDαt x(t) =

1Γ(−α)

∫ t

a(t− τ)−α−1x(t)dt, (1.14)

L’espressione ottenuta ricorda molto da vicino una convoluzione, in- . . . poi trasformando,sia con laconvoluzione diRiemann-Liuoville. . .

fatti basta definire il fattore dei pesi come semplice funzione esponen-ziale, ad esempio scrivendo f (t) ≡ t−α−1, per poter riscrivere

aDαt x(t) =

1Γ(−α)

f (t)⊗ x(t),

similmente al procedimento illustrato nel § 1.4.2 nella pagina 13.La trasformata gode della proprietà di equivalenza tra convoluzione . . . sia con la

proprietà dellaconvoluzione dellatrasformata. . .

e prodotto (1.24), si può quindi riscrivere la (1.14) come prodotto:

L {aDαt x(t)} =

1Γ(−α)

L{

t−1−α}

L {x(t)} = sαL {x(t)}

per ottenere il 2° membro dell’equazione mostrata; l’uguaglianza tra il1° ed il 3° membro invece, discende dalla proprietà (1.26) di differinte-grazione della trasformata.

Se infine si ricorda che la trasformata dell’impulso è unitaria, bastasostituire L {x(t)} = δ(t) ⇒= L {δ(t)} = 1, per ottenere dall’ugua- . . . l’uguaglianza

delle soluzioni. . .glianza tra il 2° ed il 3° membro la relazione:

L{

t−1−α}

Γ(−α)= sα,

da cui antitrasformando rispetto al tempo (la funzione gamma non . . . permette ditrovare la rispostaimpulsiva dellafrattanza!

interviene poiché dipende solo da α), si ottiene

L −1 {sα} =t−1−α

Γ(−α)u(t) (2.6)


che rappresenta la risposta impulsiva di un sistema sα cercata, ed incui il gradino u(t) si aggiunge per ricordare che l’integrale di Laplaceè tra 0 ed ∞.

Derivatore e rappresentazioni grafiche

Come accennato nell’introduzione, il termine monomio corrispondead un capacitore frazionario, cioè ad una frattanza. Si intuisce come, ana-logamente al caso continuo, un termine del genere si possa sfruttareper calcolare un differintegrale.

Consideriamo un sistema costituito proprio da un termine (2.5); po-nendo in ingresso un segnale x(t), mi aspetto che l’uscita y(t) sia parialla derivata dell’ingresso. È possibile calcolare l’uscita come di con-sueto, con il prodotto in frequenza o la convoluzione nel tempo (1.24):

y(t) = L −1 {sα} ⊗ x(t) Y(s) = sαX(s),

basta poi ricordare la (1.26) per capire come sia possibile calcolare laderivata frazionaria di una funzione, come convoluzione tra la fun-zione e la risposta impulsiva (2.6) del sistema sα; in realtà si ottienein questo modo, come era ovvio, la definizione di differintegrale diRiemann Liouville (1.14).

Per operare su sistemi di ordine non intero, non è possibile sfruttarei comandi del controlsystem toolbox. Ad esempio tentando di esegui-re il comando bode su un sistema di ordine frazionario, si ottiene ilmessaggio di errore:

1 ??? Error using ==> lti.mpower at 34

2 Power must be an integer.

Sono stati quindi sviluppati degli appositi toolbox per Matlab, in par-ticolare sono noti Ii toolbox crone e ninteger; quest’ultimo tra l’altro ègratuito1, sebbene sia meno sviluppato, e dispone di un breve articolodi documentazione scritto dal suo autore [42].

2.2.2 Termine binomio

Un’altra funzione notevole per lo studio dei sistemi è la seguente:

H(s) =1

sα − aper α > 0. (2.7)

I poli dell’equazione si possono ottenere dall’equazione sα = a, chepresenta infinite soluzioni sulla circonferenza di raggio α

√|a|. In gene-Solo se 0 < α ≤ 1 i

poli sono sullasuperficie principale

di Riemann!

rale non è possibile in assicurare l’esistenza di un polo sulla superficieprincipale di Riemann, è possibile farlo invece se si assume 0 < α ≤ 1.

Risposta impulsiva

Il procedimento per il calcolo della risposta impulsiva è stato già tratta-to nel § 1.6 nella pagina 18, anzi, la (1.42) permette di calcolare la rispo-sta impulsiva anche per termini molto più complessi rispetto alla (2.7),

1 Il toolbox è disponibile sul sito di scambio Matlab, all’indirizzo http://www.mathworks.com/matlabcentral/fileexchange/8312




per la quale basta la soluzione semplificata (1.43). Per l’inversionedunque è sufficiente la

L −1{

1sα ∓ a

}= tα−1Eα,α(±tαa), (1.43)

che fornisce una risposta impulsiva, che per α intero si riduce allarisposta impulsiva dei sistemi interi, con la funzione di Mittag-Lefflerche risulta pari all’esponenziale, come è ovvio che sia!

2.2.3 Funzioni più complesse

Sistemi più complessi, con un certo numero di poli e zeri non tuttinulli, si possono ricondurre tipicamente ai casi più semplici, §2.2.1 e§2.2.2. Cerchiamo di capire come.

Analizzando la (2.4) è possibile distinguere i casi: Esprimo gliesponenti comemultipli interi di unsolo esponentefrazionario.

1. Gli esponenti αn sono numeri razionali. In tal caso, dopo aver tro-vato il minimo comune multiplo α, indichiamo tutti gli esponenticome suoi multipli, nα per n = 0, 1, 2, . . . , N;

2. Gli esponenti α sono numeri reali, ma multipli di un valore dibase α ∈ (0, 1). Anche in questo caso è possibile scrivere gliesponenti come multipli di α, quindi ancora una volta αn = nαper n = 0, 1, 2, . . . , N.

Nel primo caso esiste un procedimento semplificato per ricavare larisposta impulsiva nel dominio del tempo, come mostrato da Ortiguei-ra ad esempio in [28, 30], nel secondo caso occorre invece la trasforma-ta della funzione di Mittag-Leffler [34], oppure la (2.6). In questo testoconsideriamo solo la seconda ipotesi, i due casi sono quindi riuniti, laprocedura da adottare è la stessa per entrambi ed è la seguente:

1. Trasformiamo H(s) in H(w), sostituendo sα con w. Se H(w) nonfosse una frazione propria, sarebbe necessario estrarre la parteimpropria come somma polinomiale;

2. Il polinomio a denominatore della H(w) è detto polinomio indi-ciale. Espandere H(w) in frazioni parziali, utilizzando i soli zeridel polinomio indiciale che giacciono sulla superficie di Riemannprincipale.

A questo punto non resta che invertire termine a termine per avere larisposta impulsiva del sistema!ESEMPIO 1. Risoluzione di un semplice sistema frazionario. Cerchiamo lasoluzione di un sistema descritto dalla seguente equazione differenziale (daun esempio del [34]) a coefficienti non interi, con le condizioni iniziali date:

{D1/2 f (t) + a f (t) = 0 per t > 0,d−1/2 f (t)

∣∣t=0= C

(2.8)

Applicando la (1.26) alla prima equazione del sistema, termine a termine,otteniamo:

L{

D12 f (t)

}= s

12 F(s)−

0

∑k=0

skD12−k−1 f (t)

∣∣t=0 = s

12 F(s)− C,


per il primo termine, in cui si è fatto uso della condizione iniziale impostadalla (2.8), e

L {a f (t)} = aF(s),

per il secondo. Rimettendo insieme i pezzi, si ha:

L{

D12 f (t) + a f (t)

}= s

12 F(s) + aF(s)− C = F(s)

(s

12 + a

)− C = 0,

dalla quale ricaviamo infine

F(s) =C

s12 + a

.

Il procedimento è applicabile in generale ad una qualsiasi equazione diffe-renziale, dalla quale si otterrà sempre un’equazione con potenze ordinate dis di termini F(s), dalla quale esplicitare infine una funzione razionale. Lapresenza di esponenti frazionari non cambia le considerazioni sulla formadelle funzioni associate al sistema svolte nel caso di esponenti interi!

In questo caso ho ottenuto una funzione elementare invertibile conla (1.42), per funzioni più complesse sarebbe stata necessaria una scomposi-zione in frazioni parziali, in questo caso invece basta riconoscere nella (1.42)che:

• l’esponente (k + 1) è pari ad 1, quindi k = 0 e l’ordine di derivazionedi M-L è 0;

• α = 1/2;

• l’esponente α− β è deve essere 0, quindi α = β;

per ottenere in fine la formula cercata:

f (t) = Ct−12 E 1

2 , 12(−at

12 )

2.3 STABILITÀCosì come per i sistemi interi, si considera la stabilità dei sistemi noninteri [28] in senso Bounded Imput Bounded Output (BIBO), cioè ilcircuito si dice stabile, se ad ingressi limitati corrispondono uscitelimitate.

Consideriamo un semplice sistema che rappresenti, ad esempio, lafunzione di trasferimento di un condensatore frazionario, sα, con sappartenente alla superficie principale di Riemann. In base al valoredell’esponente, si avranno differenti comportamenti del sistema:

• per α > 0 il sistema sarà chiaramente instabile, con un diagram-ma di Bode che cresce con pendenza di 20α dB/dec;

• per −1 ≤ α < 0 la risposta impulsiva si mantiene sempre limi-tata nel tempo, il sistema sarà stabile in senso stretto, così co-me risulta essere stabile il sistema ordinario che si ottiene perα = −1;

• per α < −1 otteniamo nuovamente un sistema instabile, la cuirisposta impulsiva per t→ ∞, tende ad infinito!


Si comprende quindi il motivo per il quale, anche un sistema polino-miale, dovrà avere il termine di grado massimo αn = nα (in cui n = N),limitato, cioè Nα < 1.

Per α = 1/n, ricordando le (1.33), si osserva come il piano complesso(0 < ϕ < 2π) si trasformi nel settore (0 < ϑ < 2π/n), mentre ilsemipiano sinistro, si trasforma nel settore π/2n < ϕ < π/2n + π/n.Lo stesso avviene per i settori traslati di 2π/n · i, per cui il sistema èstabile se i poli sono confinati in uno dei settori che si mappano nelsemipiano sinistro!

Se α non è razionale, le considerazioni sono le stesse, anche in talcaso le regioni di stabilità sono quelle che si mappano nel semipianonegativo, definite dalle fasi:

π

2α + 2παi < ϕ <

π

2α + 2παi + απ, per i = 0, 1, . . . , N − 1.

2.4 APPROSSIMAZIONII sistemi di ordine non intero permettono di modellare con più accu- I modelli frazionari

possono essere piùaccurati. . .

ratezza dei sistemi che per loro natura non sono descritti da derivateo integrali ordinari, ma frazionari. Le prime applicazioni della teoriadei sistemi di ordine non intero infatti, riguardano lo studio di sistemifisici e chimici che non possono essere ricondotti ad equazioni differen-ziali ordinarie (Ordinary Differential Equation, ODE), se non a scapitodi una consistente perdita di accuratezza del modello. Tra questi sihanno i processi elettrochimici di diffusione, i moti browniani, la vi-scoelasticità oppure il rumore di tipo 1/ f , tanto per citare i più noti,da tempo studiati [43].

Se il sistema da controllare è per sua natura frazionario, un control-lore frazionario avrà certamente migliori prestazioni, per questo sonostati sviluppati dei metodi per il dimensionamento e la progettazio-ne di vari tipi di controllori frazionari, come mostrato nella sezioneseguente, che permettono di ottenere una funzione di trasferimentofrazionaria.

Partiamo quindi da una rappresentazione di un sistema elettronicodata come funzione di trasferimento del tipo 2.4 nella pagina 21, ecerchiamo un metodo di sintesi, concentrandoci in questa parte sullarealizzazione di sistemi analogici, mentre il caso discreto sarà trattatonel prossimo capitolo.

In letteratura la strada tipicamente percorsa è stata quella dell’ap- . . . ma spesso sirealizzano conelementiconvenzionali,approssimando alcaso intero

prossimazione della funzione di trasferimento ottenuta, con una diordine intero da sintetizzare poi nel modo usuale, essenzialmente conil teorema di Darlington [4], e con elementi concentrati convenzionali,cioè, nel caso analogico, con resistenze capacità ed induttanze. La pri-ma sintesi di una rete RC con funzione di trasferimento s1/2 risale al1975, successivamente sono state studiate reti RC ed RL per la sintesidi funzioni di trasferimento con α compreso nell’intervallo (0, 1) [31].

I metodi per approssimare una funzione con esponenti frazionariin una serie di termini con esponenti interi sono molteplici, per unavisione d’insieme del problema, l’articolo [43] è un buon punto di par-tenza, per l’approfondimento di metodi specifici [9, 17, 35, 5, 13] soloper citarne alcuni!

La maggior parte dei metodi si basa sullo sviluppo in frazioni conti- approssimazionitempo-continue


nue (Continued Fraction Expansion, CFE), spesso combinato con altretecniche di pre-distorsione o condizionamento dei segnali. Le princi-pali approssimazioni nel dominio tempo-continuo:

1. Espansione in frazioni continue (CFE) [44];

2. Carlson (Newton) [9];

3. Matsuda (CFE con interpolazione) [43];

4. Oustaloup;

5. Chareff [11];

6. Minimi quadrati [43].

Le frazioni continue, rispetto ad un’espansione in serie, presentanole approssimazioniCFE sono più

utilizzate, perchépresentano una serie

di vantaggi

vari vantaggi. Infatti, le funzioni che si ottengono da un’espansioneContinued Fraction Expansion (CFE):

• convergono in genere più rapidamente delle espansioni in seriedi potenze, alla funzione che approssimano [43];

• convergono tipicamente in un dominio più ampio del piano com-plesso [43, 7];

• permettono di modellare efficacemente funzioni con i poli [43];

• nel caso discreto, originano filtri di tipo Infinite Impulse Response(IIR), che convergono più rapidamente a parità di complessitàcomputazionale, delle approssimazioni Finite Impulse Response(FIR) [5].

Vediamo quindi di che si tratta, e come l’espansione CFE possa esseresfrutta convenientemente per modellare una funzione di trasferimento.

2.4.1 Espansione in frazioni continue

Le basi della teoria delle frazioni continue sono state poste da ThomasJoannes Stieltjes (1856-1894) sul finire dell’ottocento, nel 1948 Wall[44], nel suo libro sulle frazioni continue [44], ha raccolto ed organiz-zato in maniera coerente la materia; un approccio meno teorico edorientato alle applicazioni pratiche è mostrato in [21].

Le frazioni continue possono assumere varie forme, tutte riconduci-bili alla seguente espressione2:Espressione estesa

f ≡ b0 +a1

b1 +a2

b2 +a3

b3 + · · ·

,

spesso resa con la forma più compattaforma compatta

f ≡ b0 +a1

b1+a2

b2+. . . ,

dovuta a Rogers .


È possibile esprimere in questo modo un qualsiasi numero reale, nel sviluppo CFE per inumeri realisenso che una serie infinita di frazioni continue con opportuni coeffi-

cienti a, b, converge al numero dato. Tale espansione è equivalente allarappresentazione decimale, nel senso che è unica, ma presenta alcunivantaggi. Ad esempio il numero razionale 1/3 = 0.333 . . . ha una rap-presentazione decimale infinita, in frazioni continue è data invece daitermini a1 = 1, b1 = 3, così come in generale ogni numero razionale hauna rappresentazione finita.

Le frazioni continue sono generate da un’infinita sequenza di tra- definizione rigorosadi CFEsformazioni lineari dalla variabile x alla variabile y [44], secondo le

relazioni

y0(x) = b0 + x

yk(x) =ak

bk + x, k = 1, 2, 3, . . .

di ricorrenza.Le prime n trasformazioni y0y1 . . . yn(0) = y0y1 . . . yn+1(∞), genera-

no l’n-esimo approssimante

y0y1 . . . yn(0) = b0 +a1

b1 +a2

b2 + . . .an

bn + x

che è possibile riscrivere, con l’ausilio delle relazioni di ricorrenza formule frazionariericorrenti

A−1 = 1, B−1 = 0,

A0 = b0, B0 = 1,

Ak+1 = bk+1 Ak + ak+1 Ak−1,

Bk+1 = bk+1Bk + ak+1Bk−1,

in forma di frazione semplice di grado n:

y0y1 . . . yn(x) =An−1x + An

Bn−1x + Bn,

come si dimostra facilmente per induzione [44, cap. 1].Il problema per l’espansione CFE di una funzione, consiste nella de-

terminazione dei coefficienti a, b che portano la frazione continua aconvergere alla funzione da approssimare. Sono stati sviluppati de-gli algoritmi iterativi che permettono di identificare i coefficienti peralcune classi di funzioni.

Noti i coefficienti dello sviluppo CFE, può essere comodo a volterisolvere e riportare la frazione continua in una frazione semplice, ilseguente script si occupa proprio di questo compito.

1 A(1)=1;

2 B(1)=0;

3 A(2)=b(1);

4 B(2)=1;

5

6 for n=3:7,

2 Il termine ak è detto numeratore parziale k-esimo, il termine bk è detto denominatoreparziale k-esimo, infine il loro rapporto ak/bk è detto k-esimo quoziente parziale.


7 A(n)=b(n).*A(n-1) + a(n).* A(n-2);

8 B(n)=b(n).*B(n-1) + a(n).*B(n-2);

9 end

10

11 c=[A;B]

Funzioni razionali

Consideriamo una frazione espressa dal rapporto di due polinomi,come mostrato in [21, capitolo I, §3], del tipo

f (x) =c10 + c11x + c12x2 + c13x3 + . . .c00 + c01x + c02x2 + c03x3 + . . .

.

È possibile ottenere una scomposizione CFE

f (x) = b0 +a1

b1 +a2x

b2 +a3x

b3 + · · ·

,

con un algoritmo ricorsivo, sintetizzato nelle seguenti formule:Algoritmo ricorsivoper lo sviluppo in

frazioni parziali diun rapporto di

polinomi

c20 = c10c01 − c00c11

c21 = c10c02 − c00c12

c22 = c10c03 − c00c13

· · ·c30 = c20c11 − c10c21

c31 = c20c12 − c10c22

c32 = c20c13 − c10c23

· · ·cmn = cm−1, 0 cm−2, n+1 − cm−2, 0 cm−1, n+1,

che ricordano il criterio di Routh, si osservi infatti, come ciascun ter-mine sia dato dal determinante della matrice data dall’incrocio delledue righe precedenti con la prima colonna e quella seguente al terminecalcolato.

Dalla prima colonna della tabella così ottenuta, posto il primo coef-ficiente b0 = 0, le formule

ai = ci,0 per i = 1, 2, 3, . . .

bi = ci−1,0, per i = 1, 2, 3, . . . ,

identificano la serie completa di coefficienti della CFE.Nel toolbox ninteger per Matlab, è disponibile il comando contfracf

che implementa proprio l’algoritmo descritto, come mostrato nell’ap-pendice A.

Potenze non intere

Un termine molto usato nella teoria dei sistemi frazionari è un terminedella forma

f (x) = (1 + x)α. (2.9)


Il suo sviluppo è un caso particolare delle serie ipergeometriche,

F(a, b, c; z) = 1 +abc

z +a(a + 1)b(b + 1)

2!c(c + 1)z2 + · · ·+ (a)n(b)n

(n− 1)!(c)nzn + · · ·

il cui sviluppo CFE è stato trovato da Carl Friedrich Gauss (1777 – lo sviluppo CFE delleserie ipergeometricheè stato ricavato daGauss

1855) matematico, astronomo e fisico tedesco. I simboli (a)n indicanoil fattoriale crescente, cioé (a)n = a(a + 1)(a + 2) · · · (a + n− 1), e sonodovuti a Pochammer .

I coefficienti dello sviluppo, nel caso particolare in cui F(−α, 1, 1;−z) =(1 + z)α, sono dati da

ak ={−α,

1(1 + α)1 · 2 ,

1(1− α)2 · 3 ,

2(2 + α)3 · 4 ,

2(2− α)4 · 5 , . . .

},

si noti l’alternanza tra termini pari e dispari. I coefficienti b dellosviluppo sono tutti unitari, tranne il primo che è nullo.

Nel toolbox ninteger non è presente alcuna funzione per calcolare implementazioneMatlablo sviluppo in frazioni continue della funzione (2.9), per questo è stata

scritta la funzione cfracpwr che prende in ingresso l’esponente α, ed ilnumero di termini dell’approssimazione n, per restituire il vettore deicoefficienti a dello sviluppo CFE.

1 a = cfracpwr (k,n);

2

3 i=1:(ceil(n/2)-1);

4 a(1)=1;

5 a(2)=-k;

6

7 a(2*i+1)= (i.*(i+k))./((2*i).*((2*i)-1));

8 a(2*i+2)= (i.*(i-k))./((2*i).*((2*i)+1));

Vengono calcolati per prima cosa i primi due coefficienti, poi i ter-mini dispari, infine i termini pari dello sviluppo.

Considerazioni

Dal punto di vista matematico, le funzioni di trasferimento come la 2.4nella pagina 21, sono razionali se gli esponenti sono tutti interi, altri-menti, se compaiono esponenti non interi, sono irrazionali, analoga- Il problema della

rappresentazioneillimitata

mente a quanto avviene nel passaggio dall’insieme dall’insieme deinumeri razionali Q all’insieme dei numeri reali R. Si può dire dunqueche la rappresentazione è illimitata in s, e che il caso con esponentiinteri sia solo un caso particolare con rappresentazione finita [43].

L’espansione CFE per le funzioni presenta degli ulteriori vantaggirispetto all’espansione in serie di potenze o altri metodi, sia per circuitianalogici che digitali.

Le funzioni che presentano poli sono difficilmente approssimate cor- funzioni con polirettamente dall’espansione in serie [43]: una serie di potenze non puòtendere ad infinito in un punto nel cui intorno è limitata. Le CFE

permettono invece di modellare i poli in modo naturale, perché sonopropri delle frazioni di cui fa uso lo sviluppo.

Le trasformazioni che definiscono le funzioni continue, creano una spazi trasformaticorrispondenza tra 2 variabili di due distinti spazi. Passando dal domi-nio continuo di Laplace s al dominio discreto z con una trasformazione,come la bilineare, si ha la corrispondenza tra semipiani di s e regioni


circolari di z. Le frazioni continue si prestano naturalmente a tale tra-sformazione, per definizione: infatti le trasformazioni che definisconole frazioni continue stabiliscono una corrispondenza tra due dominiappartenenti a differenti spazi [44, capitolo XV]!

Le approssimazioni ottenute presentano in generale poli e zeri in-poli e zeriinterlacciati terlacciati sull’asse reale negativo. Aumentando il grado dell’approssi-

mazione, i poli e gli zeri si avvicinano sempre di più, ed al limite, unnumero infinito di punti interlacciati si può assimilare ad una linea ditaglio della funzione [17], cioè se la funzione sα con 0 < α < 1 vienevista come un operatore, tale operatore sarà discontinuo sul semiassereale negativo del piano complesso, altrimenti è priva di poli e zeri senon quello nell’origine. La presenza di molteplici poli e zeri interlac-ciati è quindi modellata in modo molto più semplice da un sistemafrazionario.

2.4.2 Metodo di Carlson

Il metodo proposto da Carlson [9] è basato prevalentemente sul meto-do di Newton, per la ricerca approssimata della radice di una funzionetramite un algoritmo iterativo basato sulle derivate.

Ricordiamo in breve il metodo di Newton: supponiamo sia x la radi-metodo di Newtonce quadrata di un numero a, possiamo scrivere l’identità x2 = a, in cuix, la radice cercata è incognita. Trattando il problema come la ricercadegli zeri di una funzione, definisco f (x) ≡ x2− a, l’obiettivo a questopunto è trovare il valore di x per il quale f (x) = 0. Partendo da unpunto arbitrario3 x0 mi sposto verso la radice lungo la tangente del pun-to dato, cioè cercando lo zero della tangente invece che della funzione,e ripartendo poi da ciascun punto raggiunto reitero l’algoritmo.

Sfruttiamo quindi la ricorrenza

xi = xi−1 −f (x)f ′(x)

∣∣∣∣x=xi

per trovare l’approssimazione H(s) di una funzione [G(s)]α con espo-nente non intero, in modo del tutto analogo. Tipicamente si assumecome valore iniziale H0(s) = 1, poi si procede cercando gli zeri dellafunzione f = H1/α − G.

Riepilogando quindi, si procede nel seguente modo:

[G(s)]α equazione nota

f (x) = x1α − G(s) = 0 funzione ausiliaria

H0(s) = 1 condizioni iniziali

Hi(s) = Hi−1 −ff ′

∣∣∣∣xi=Hi−1

algoritmo di ricorsione

Per ottenere una maggiore rapidità di convergenza è possibile pre-tecniche dipre-distorsione:

Aggiungere unapiccola illustrazione

distorcere la funzione ausiliaria f (x) in modo opportuno. Ad esempio,se 1/α ≡ q è dispari, posso considerare la funzione ausiliaria fm(x) =(xq − G)/xm con m =, che possiede le stesse radici di f , ma è piùlineare intorno alla radice.

3 Qui il metodo viene vagamente ricordato per poterne emulare il principio, le condizionidi convergenza esulano dallo scopo per cui non sono trattate.


Intuitivamente più è elevato il grado q maggiore è la curvatura del-la funzione, che quindi convergerà più lentamente! Analiticamente, lalinearità nell’intorno della radice si ha con l’annullamento della deriva-ta seconda della funzione, condizione dalla quale si ottiene un valoreconveniente dell’esponente m.ESEMPIO 2. Approssimazione di una capacità frazionaria. Supponiamo che ilsistema da approssimare G(s) sia

G(s) =(

1s

) 13

,

cioè un capacità di ordine 1/3, per cui q = 3.Per una maggiore rapidità di convergenza applichiamo una distorsione,

con m = (3− 1)/2 = 1, si ottiene la funzione

fm(x) = x2 − 1sx

,

e la sua derivata sarà

f ′m(x) = 2x +1

sx2 .

Dalla condizione iniziale H0 = 1 l’algoritmo si sviluppa

H0 = 1

H1 = H0 −fm

f ′m

∣∣∣∣x=H0

= 1− x2 − (sx)−1

2x + (sx2)−1

∣∣∣∣x=1

=s + 2

2s + 1.

Le iterazioni successive crescono rapidamente di grado. Il quadrato deltermine H1(s) sarà infatti

H1(s)2 =s2 + 4s + 4

4s2 + 4s + 1,

che sostituito nella formula di ricorsione ci da la seguente:

H2 = H1 −x2 − (sx)−1

2x + (sx2)−1

∣∣∣∣x=H1

=s5 + 24s4 + 80s3 + 92s2 + 42s + 4

4s5 + 42s4 + 92s3 + 80s2 + 24s + 1

approssimazione di 5° grado in due sole iterazioni.

2.4.3 Metodo di Chareff

Il Metodo descritto da Chareff [11], introduce un polo o uno zero,ogni volta che la funzione approssimata si scosta del limite massi-mo ammesso dalla funzione da approssimare. La pendenza arbitra-ria arbitraria −20α dB/decade di un polo di ordine α, si assimila aduna linea spezzata con due sole possibili pendenze, −20 dB/decade e0 dB/decade.

Aumentando il numero di poli e di zeri nell’intervallo di approssima-zione, aumenta l’accuratezza raggiunta cioè si riduce lo scostamentoy dalla funzione. Al di fuori dell’intervallo in cui è presente l’appros-simazione la pendenza della curva può essere un multiplo intero di20 dB/decade, o al limite essere nulla.

Considerando la funzione passabasso

G(s) =1(

1 + sp

)α , (2.10)


il cui andamento, quasi esatto in bassa frequenza, è dato dalla funzioneapprossimante

H(s) =∏N−1

i=0

(1 + s

zi

)∏N

i=0

(1 + s

pi

) (2.11)

che ha il primo polo proprio in prossimità di p p, mentre in alta fre-quenza, esternamente all’intervallo approssimato, cioè dopo l’ultimopolo o zero, si può scegliere se approssimare con pendenza piatta odi −20 dB/decade, sebbene il secondo caso sia nella maggioranza deicasi più conveniente.

Considerando come criterio, la limitazione dell’errore al di sotto diIl metodo di Chareffsi basa sulla

minimizzazione delladistanza tra la

funzioneapprossimata e la

funzioneapprossimante

un valore y scelto, è possibile posizionare ciascun polo e zero, con lesemplici relazioni:

p0 = p10y

20α , z0 = p010y

10(1−α) , (2.12)

P1 = z010y

10α , z1 = p110y

10(1−α) ,...

...

pN = zN−110y

10α , zN−1 = 10y

10α .

Osservando poi che i termini sono in rapporto geometrico, cioè ilrapporto tra ciascun termine ed il precedente è costante, si può scrivere

pi = (ab)i p0, per i = 1, 2, 3, . . . , N − 1,

zi = (ab)iap0 per i = 1, 2, 3, . . . , N − 1,

(2.13a)

(2.13b)

in cui a, b, si ottengono dai rapporti

a =z0

p0=

z1

p1= · · · = zN−1

pN−1= 10

y10(1−α) , (2.14a)

b =p1

z0=

p2

z1= · · · = pN

zN−1= 10

y10α , (2.14b)

z1

z0=

z2

z1= · · · = zN−1

zN−2= ab,

p1

p0=

p2

p1= · · · = pN

pN−1= ab,

che danno tutti luogo a serie geometriche.Il numero di punti necessario per approssimare un termine semplice

del tipo (2.10) si ottiene considerando che la spaziatura dei poli e deglizeri è in progressione logaritmica, per cui considerando la massimafrequenza ωMax superiore da raggiungere, questa cadrà nell’intervallo

(N − 1) log(ab) < log(

ωMax

p0

)< N log(ab), (2.15)

per cui l’estremo superiore sarà il numero minimo di punti che garan-tisce un’approssimazione sufficiente.


ESEMPIO 3. Approssimazione di un termine frazionario semplice. Consideria-mo un termine frazionario semplice del tipo mostrato in formula (2.10), perp = 1 e α = 1/3, cioè

G(s) =1

(1 + s)1/3 ,

e cerchiamone un’approssimazione con il metodo di Chareff che non sia maipeggiore di y = 3 dB fino ad una frequenza ωMax = 105 rad s−1.

Approssimando con il metodo di Chareff, per prima cosa bisogna calcolarei rapporti a, b, con le (2.14) si ottiene a = 2.8184, b = 7.9433. Dal loro prodottoa · b = 22.3872 si ottiene il termine base per calcolare poli e zeri, che sarannoin numero

N =

log(

ωMaxp0

)log(ab)

= 4,

come si ottiene dall’inversione della seconda disuguaglianza della (2.15),approssimata all’intero superiore (come indicato dalle parentesi de).

Il calcolo delle (2.13) si implementa molto semplicemente in Matlab conle istruzioni

1 i=0:N;

2 p=ab.^i * p0;

3 i(N)=[];

4 z=ab.^i * a * p0;

in cui p0 è il primo polo calcolato con la formula esplicita (2.12). Dal-l’algoritmo si ottengono i valori dei poli e degli zeri della formulaapprossimante (2.11), che qui riportiamo per completezza:

p 2.818 38 63.095 1412.53 31 622.7 707 945z 7.943 28 177.827 3981.07 89 125.0

2.4.4 Metodo di Oustaloup

Il metodo di Oustaloup [31] è simile al metodo di Chareff, con variazio-ni delle formula per l’approssimazione dei coefficienti, per maggioriinformazioni si rimanda al testo in bibliografia.

2.5 SISTEMI DI CONTROLLOIn letteratura sono elencati numerosi esempi di sistemi fisici che sonodescritti naturalmente da equazioni differenziali frazionarie! Progetta-re un controllore per un processo di questo tipo permette di beneficia-re di due importanti vantaggi se il controllore è un controllore linearefrazionale!

• Una descrizione più realistica del modello, che non approssimi isistemi integro-differenziali di ordine frazionario con sistemi diordine intero, permette di progettare un controllore più accuratoe più stabile, in grado di prevedere con maggiore accuratezza l’u-scita del sistema ed adattare di conseguenza l’azione di controllo,che risulterà naturalmente più robusta;

• Disporre di qualche grado di libertà in più, così come avviene


per controllori di ordine non intero, permette di ottenere unafunzione di trasferimento che approssima meglio la caratteristicadesiderata. Ciò sarà particolarmente vero, se la funzione di tra-sferimento desiderata in frequenza, non ha pendenze multiple di10 o 20 decibel per decade;

Tuttavia, i metodi per lo sviluppo di controllori di ordine non intero,non possono ancora competere con la facilità di sviluppo e l’ampioventaglio di tecniche esistenti nel caso di controllori ordinari [7], adesempio non esiste un metodo equivalente al luogo delle radici!

Consideriamo un controllore che possa scriversi nella forma [cfr. 4]:

G(s) = K +V0

sλ+ sδT(0), (2.16)

in cui usualmente i parametri λ, δ, sono interi (e unitari), che generaliz-ziamo al caso di esponenti reali. Tale struttura permette di realizzareun filtro con molta più flessibilità rispetto a quanto offerto da un’archi-tettura a coefficienti unitari o comunque interi, e permette di realizzareun controllore flessibile, che può adattarsi meglio alle dinamiche deisistemi controllati, specie se frazionari!

Nello spazio di stato un sistema si può esprimere nella forma:

g(t) =

{Dαx(t) = Ax(t) + Bu(t)y(t) = Cx(t) + D(u),

Quando α non è intero unitario, ma si ha 0 < α ≤ 1, la soluzione èesprimibile nel modo consueto, tramite una sostituzione dell’esponen-ziale con la più genereale ed equivalente funzione di Mittag-Leffler:

x(t) = Φ0(t)x0 +∫ t

0Φ(t− τ)Bu(τ)dτ

in cui

Φ0(t) = Eα(Atα),

Φ(t) =∞

∑k=0

Akt(k+1)α−1

Γ[(k + 1)α]

come è possibile dimostrare seguendo le indicazioni di [20], dalle qualiè anche possibile studiare l’osservabilità del sistema.

È possibile tracciare i diagrammi di Bode della (2.16) nel modo usua-le, con gli strumenti forniti. Dai grafici ottenuti è possibile ricavare lastabilità del sistema nel modo usuale, ad esempio con il criterio diNyquist, per funzioni implementate in grado di tracciare i diagrammicitati, si consulti l’appendice A nella pagina 83 .

2.5.1 Controllore lead-lag

Sono i classici sistemi di controllo con reti di ritardo o di anticipo, madi ordine non intero [48]. Se tipicamente un rete anticipatrice assume laforma

C =1 + τas1 + τa

mas

,


nel caso di un controllore non intero si avrà una funzione compensa-trice del tipo

C =

(1 + τas1 + τa

mas

,

)α

, (2.17)

in cui τa è la costante tempo, la “posizione” temporale in cui agiscela rete, mentre ma è il margine di guadagno, cioè l’anticipo di fasepermesso dalla rete. Ovviamente si può costruire una rete analoga cheritarda la fase.

La presenza dell’esponente α nella (2.17) influisce sia sui moduli,che avranno una pendenza massima di α20dB/dec, sia sulle fasi, chepossono avere uno sfasamento massimo di α90°.

La funzione (2.17) corrisponde ad un derivatore/integratore frazio-nario limitato in frequenza, che sta alla base dei sistemi di controlloContrôle Robuste d’Ordre Non Entier (CRONE) [26].

Il metodo risulta adeguato per piccole correzioni di sistemi che giàpossiedono una discreta stabilità, risulta invece inefficace per sistemifortemente instabili, analogamente a quanto avviene nel caso intero.

2.5.2 Controllori di tipo PID

Sono i controllori costituiti dal parallelo di tre blocchi, un blocco pro-porzionale (p) che realizza una funzione di amplificazione idealmentecostante Kp, un blocco integratore (i) che taglia le alte frequenze 1/s, edun blocco derivatore (d) che taglia le basse frequenze s [10, 48].

Un controllore Proportional Integral and Derivative (PID) assumequindi la forma

C(s) = Kp(1 + kis−1 + kds),

con le tre funzioni specificate, ciascuna col suo guadagno, ed i terminiin s tutti di ordine intero.

Nel caso frazionario la generalizzazione per i sistemi di tipo ProportionalIntegral and Derivative (PIλDµ) è data dalla semplice estensione

C(s) = Kp(1 + kis−λ + kdsµ), (2.18)

in cui ciascun termine in s assume un esponente frazionario.Come evidenziato da [10], l’ostacolo maggiore nella realizzazione

di controllori PIλDµ risiede nelle difficoltà di tuning, che richiedonometodi soluzione non lineari alquanto complessi, specialmente con ilcrescere dei parametri dovuto all’introduzione degli esponenti frazio-nari λ e µ. Lo sforzo realizzativo è giustificato dalla grande flessibilitàofferta in fase di design del controllore.

2.5.3 Controllori di tipo TID frazionari

Sono una variante dei controllori PID che offre maggiore versatilità, incui il blocco proporzionale è sostituito da una funzione rimodellante,cioè da un blocco di tipo tilted della forma 1/s(1/n).

La funzione di trasferimento tipica del metodo è quindi

C(s) = Kp

(1

s(1/n)+ kis−λ + kdsµ

),

con i vantaggi e le complicazioni che il metodo comporta.


2.5.4 Controllori CRONE

Tecnica di controllo sviluppata dal team capitanato dal ricercatore fran-cese Oustaloup. Attualmente i controllori Contrôle Robuste d’OrdreNon Entier (CRONE), sono tra i più robusti, e sono disponibili di tregenerazioni. Per la progettazione è stato sviluppato un apposito toolper Matlab a pagamento, che prende appunto il nome di crone. Permaggiori approfondimenti sul metodo si rimanda a [33].

3 S I S T E M I T E M P O - D I S C R E T I

Lo sviluppo di sistemi frazionari discreti non è attualmente avanza-to quanto la rispettiva controparte continua: le applicazioni sono at-tualmente poco numerose, quelle presenti riguardano principalmentel’implementazione digitale di sistemi di controllo, tuttavia i risultatiottenuti sono piuttosto promettenti, tanto da incoraggiare lo sviluppoe la ricerca!

La caratteristica principale dei sistemi digitali non interi, detti a ritar-do frazionario, è costituita dalla presenza di ritardo o anticipo (lead) noninteri, ma frazionari. Rispetto ai sistemi discreti ordinari, richiedonol’interpolazione dei segnali, e trovano svariate possibilità applicativenell’ambito dell’elaborazione della voce, della musica, di array e moltialtri.

I sistemi digitali si possono ottenere come sistemi a ritardo frazionarioa partire da un’approccio completamente digitale del problema, oppuresi possono ottenere dai sistemi frazionari tempo-continui, per discre-tizzazione dei differintegrali o della variabile s di Laplace. In tal caso,un solo differintegrale da origine ad un sistema che richiede molteplicicampioni del segnale in elaborazione.

La formulazione matematica che si ottiene risulta analoga alla for-mulazione di sistemi di ordine intero. Come avviene nel caso tempo-continuo però, si ha a che fare con rappresentazioni irrazionali, oppureinfinite, subentrano quindi molteplici difficoltà matematiche, assentinel caso intero.

3.1 INTRODUZIONEUn solo articolo proposto da Ortigueira [29], per quanto ci è stato possi-bile verificare, offre una visione d’insieme dei sistemi discreti di ordinenon intero a ritardo frazionario. Altri articoli, non molti in verità, affron-tano argomenti specifici, come sistemi di controllo discreti o tecnichespecifiche di discretizzazione ed approssimazione, che possono essereapplicate anche a sistemi discreti derivati da sistemi tempo-continuo.Per quel che riguarda le tematiche legate al ritardo o anticipo nonintero, un’analisi dettagliata è stata proposta da Laakso et al. [23].

Cerchiamo di rivedere ed espandere le tematiche legate ai sitemi aritardo non intero, come se fossero sistemi distinti e definiti, piuttostoche studiarli nel modo usuale, basato essenzialmente sul ricampiona-mento, sebbene questo sia sempre uno strumento fondamentale.

Nei sistemi discreti il dominio del tempo è rappresentato dall’insie- ritardo/anticipofrazionariome dei numeri interi n, con quali sono ordinati i campioni del segnale.

I campioni del segnale si succedono ad intervalli di tempo regolari, lacui durata è indicata comunemente con T. In alcuni casi può essere ne-cessario conoscere il valore assunto tra due punti di campionamentoattigui, in altri può essere necessario sincronizzare due segnali cam-pionati con lo stesso periodo ma in istanti diversi. Quando il ritardo

39

40 SISTEMI TEMPO-DISCRETI

o l’anticipo possono essere anche diversi dal periodo, allora ho unsistema discreto di ordine non intero.

Le generalizzazioni non sono così immediate, il concetto di siste-Il concetto di sistemadiscreto non

intero. . .ma discreto non intero sfugge ad una interpretazione chiara, non èneppure ben chiaro come definire un sistema di questo tipo. Possiamosupporre in prima analisi di staccarci dall’equazione differenziale (2.2),e considerare il suo equivalente discreto di ordine intero:

N

∑k=0

ak y[n− k] =M

∑k=0

bk x[n− k], (3.1)

e supporre che il suo equivalente di ordine non intero, sia della forma

N

∑k=0

aky[n− αk] =M

∑k=0

bkx[n− αk]. (3.2)

Una tale equazione sarebbe giustificata da una possibile corrispon-. . . si può basaresull’equivalenza tra

termini in sα e zα. . .denza tra termini in s della forma sα e termini discreti della forma zα.Se termini del tipo zα sono distintivi dei sistemi di ordine non inte-ro a ritardo frazionale, allora ci si aspetterebbe di operare su sistemidescritti da equazioni del tipo

N

∑n=0

anzαn Y(z) =M

∑m=0

bmzαm X(z),

e l’equazione (3.2) si potrebbe ottenere antitrasformando proprio l’e-quazione in z mostrata.

Cominciamo quindi, proprio dallo studio dei termini zα, ma primauna breve considerazione. Si può osservare come un’integrazione diordine non intero possa essere posta in relazione con una integrazioneordinaria di una funzione, in cui la variabile di integrazione è scala-. . . oppure può

sorgere dal concettodi scaling implicitonelle definizioni dei

differintegrali

ta con continuità [16, pagina 82]. Brevemente, i passaggi indicati siapplicano a partire dalla definizione di Riemann-Liouville (R-L) (1.14),definendo una funzione ausiliaria di scaling

g(y) ≡ xα − (x− y)α

Γ(α + 1), (3.3)

che sostituisce la variabile ausiliaria y, cioè y→ g(y), per cui si ottiene

aD−αx =

∫ x

af (y)dg(y),

nella forma di integrale ordinario!A partire da questa idea pensiamo che sia possibile ricavare un’e-

quazione alle differenze rappresentativa di un sistema di ordine nonintero. Si potrebbe campionare la funzione considerata ad intervalliy regolari, cioè con un campionamento variabile secondo la funzionecampionamento

variabile g(y) (3.3) e trattare il sistema come un sistema ordinario, oppure cam-pionare regolarmente la funzione in esame e poi trovare le differenze didifferenze non intereordine non intero che permettano di passare ad un sistema discreto diordine frazionario.

3.2 RITARDO ELEMENTARE ZA

Cerchiamo di capire quindi per prima cosa, che significato potreb-Impostiamo l’analisi,a partire

dall’antitrasformatadi un termine zα

3.2 RITARDO ELEMENTARE ZA 41

be avere un termine zα, che in base alle considerazioni di carattereintuitivo proposte, si potrebbe assimilare ad un ritardo frazionario.

Ritardo intero

Ricordiamo alcune note relazioni di ordine intero, che cercheremo digeneralizzare per ordine qualsiasi. Un ritardo intero nel tempo, corri-sponde in frequenza, per una ben nota proprietà della trasformata Z 1,al prodotto della trasformata del segnale per un termine un termine Un termine z−n

rappresenta unritardo intero. . .

z−n. Si ha cioè l’equivalenza:

z−n0 X(z−1)⇔ x[n− n0]. (3.4)

Un termine della forma z−n0 corrisponde alla funzione di trasferimen-to di un sistema discreto elementare. Per definizione la funzione ditrasferimento H(z) si ottiene dall’uscita del sistema stesso quando iningresso è presente un impulso. L’uscita si può ottenere dal prodottoin frequenza

H(z) = Z {δ[n]} · z−n0 = 1 · z−n0 ,

piuttosto che dalla convoluzione nel tempo, ricordando che la trasfor-mata Z di un impulso è pari ad 1, quindi invertendo la funzione di . . . perciò ha una

risposta impulsivapari ad una delta diDirac traslata. . .

trasferimento

h[n] = Z −1 {H(z)} = Z −1 {z−n0}

= δ[n− n0], (3.5)

si ricava che la risposta impulsiva di un sistema z−n0 è pari ad unimpulso, traslato di n0 campioni.

Ritardo frazionario

Le cose si complicano considerando un analogo sistema che indiche-remo con X(z) = z−α con α ∈ R, che ipotizziamo sia un sistema conritardo non intero α. Se ne calcoliamo la trasformata Z inversa, con ladefinizione integrale bilatera di antitrasformata,

x[n] = Z −1 {X(z)} =1

2π j

∮γ

X(z)zn−1dz, n ∈ Z (3.6)

appoggiandoci alla teoria delle funzioni olomorfe, non è possibile ot-tenere alcuna funzione che sia causale o anticausale!

È noto dalla teoria delle funzioni olomorfe [39], che la curva chiusaγ deve essere interamente contenuta nella regione di analiticità dellafunzione, e che deve contenere al suo interno tutti i poli.

Se il sistema X(z) = z−n0 è di ordine intero, ho un solo polo inz = 0 di molteplicità n0, per cui è sufficiente che la curva γ racchiudal’origine del piano complesso C. Al contrario, se il sistema X(z) = z−α

è di ordine non intero, dovranno essere considerati molteplici fogli diRiemann sovrapposti, congiunti attraverso delle linee di taglio, o branchcut, in z = 0 ho un branch point, cioè un punto di taglio (si rilegga il § 1.5nella pagina 15 per maggior chiarezza). Ciascuna linea di taglio siscelga, comunque si avrà una semiretta che parte dall’origine e taglia

1 Il simbolo Z viene dall’iniziale del matematico ed informatico russo (Azerbaijan, 1921)Lotfali Askar-Zadeh, docente di informatica all’università di Berkley, che nel 1952

introdusse la trasformata Z .


la curva di integrazione γ, non sarà perciò possibile identificare dueregioni circolari del tipo |z| < 1, |z| > 1, e per questo non non èpossibile imporre la causalità attraverso una scelta opportuna dellaregione di taglio [29].

D’altra parte si osserva [29] come la funzione z−α per α non intero,La trasformata Z −1

di un termine zα nonesiste in senso

classico. . .

non abbia alcuno sviluppo convergente in serie di MacLaurin2 o di Lau-rent3, per cui non esiste alcuna sequenza causale o no, che abbia z−α

come trasformata Z ! Lo sviluppo in serie esiste invece se |z| = 1.

Trasformata Z a partire dalla DTFT

Semplifichiamo quindi le cose, calcolando la trasformata Z solo sulla. . . per cui calcoliamoZ −1 solo sulla

circonferenzeunitaria. . .

circonferenza unitaria, cioè occorre imporre la condizione

z⇔ ejω,

che ci ricorda il passaggio dalla trasformata Z alla Discrete TimeFourier Transform (DTFT).

Sulla circonferenza unitaria la definizione di antirasformata (3.6) sisemplifica,

x[n] =1

2π j

∫ π

−πX(ω)ejωn−1

dejω,

esplicitando poi il differenziale rispetto ad ω, si ottiene la relazionedejω = ejω jdω che ci permette di completare la sostituzione, ottenendoinfine:

x[n] =1

2π j

∫ π

−πX(ω)ejωne−jωejω jdω =

12π

∫ π

−πX(ω)ejωndω, (3.7)

che coincide proprio con la definizione di DTFT inversa, come ci aspet-. . . cioè calcoliamo laDTFT inversa. . . tavamo.

Antitrasformiamo con la (3.7) il sistema considerato, X = z−α, chein frequenza è X(jω) = e−jωα, dopo alcuni passaggi

h[n] =1

2π

∫ π

−πejω(n−α)dω =

12π

[ejω(n−α)

j(n− α)

]π

−π

=

=1

2π

[cos[ω(n− α)] + j sin[ω(n− α)]

j(n− α)

]π

−π

=

=1

2π

cos[π()] + j sin[π()]− cos[−π()]− j sin[−π()]j(n− α)

=,

proseguendo poi, grazie alle relazioni goniometriche sugli angoli op-posti (nota che sono stati omessi in parte gli argomenti (n− α) dellefunzioni sin e cos per brevità), si semplificano le funzioni cos, mentrele funzioni sin si sommano. . . otteniamo infine

un rispostaimpulsiva di tipo

sinc h[n] =sin[π(n− α)]

π(n− α), (3.8)

2 Colin MacLaurin (1698–1746), matematico scozzese. La serie del matematico ingleseBrook Taylor (1685-1731), pubblicata nel 1715 nel Methodus Incrementorum Directa etInversa, calcolata attorno all’origine, prende il nome di serie di MacLaurin.

3 Pierre Alphonse Laurent (1813–1854). La sua generalizzazione in C della serie di Taylor,che fa uso del teorema di Cauchy per il calcolo dei coefficienti, risale al 1843.


riusciamo infine ad ottenere la risposta impulsiva in frequenza del ri-tardo frazionario z−α. Quando il ritardo α assume un qualsiasi valoreintero, il seno della (3.8) si annulla, per cui l’equazione si riduce ad unadelta (3.5). Per valori non interi invece, la risposta impulsiva ha un du-rata illimitata, e non è neppure causale dal momento che si estendesia verso −∞ che verso +∞! Non è quindi implementabile in alcunaapplicazione in tempo reale. In questo modo si può trovare un’al-tra spiegazione, più intuitiva rispetto a quelle date all’inizio di questoparagrafo, che giustifichi la mancanza di una trasformata Z [23].

Ritardo frazionario ed interpolazione

È possibile rielaborare in forma diversa l’equazione trovata (3.8), ap-plicando le formule di sottrazione, diventa:

sin(πn) cos(πα)− sin(πn) cos(πα)π(n− α)

=− sin(πn) cos(πα)

π(n− α),

osservando poi che il sin(πn) assume ciclicamente i valori 0, 1 si ottie-ne:

h[n] =sin(πα)

πα

(−1)n

1− nα

, n ∈ Z, (3.9)

Possiamo calcolare l’uscita del sistema per un ingresso qualsiasi x[n], . . . quindi unricampionamento delsegnale!

come la convoluzione tra la risposta impulsiva del sistema h[n] trovatanell’equazione (3.8), e l’ingresso, ottenendo:

xn−α =+∞

∑m=−∞

x[m]sin[π(n−m− α)]

π(n−m− α)α ∈ (0, 1), n ∈ Z, (3.10)

equazione che ricorda molto da vicino il teorema del campionamen-to di Shannon-Whitaker. Abbiamo ricavato quindi, che, un ritar-do frazionario corrisponde ad una interpolazione del segnale ad unafrazione α del periodo, come ci attendevamo!

3.2.1 Differenze frazionarie

Abbiamo visto che un sistema del tipo z−α introduce un ritardo fra-zionario nella catena di elaborazione, per cui questo si comporterà co-me un’interpolatore la cui risposta impulsiva è stata ricavata ed è datadall’equazione (3.8) o dalla (3.9).

Per ottenere qualcosa che si avvicini ad un’equazione alle differenzefinite (3.1), cerchiamo di dare una definizione di derivata frazionariache sia basata sulla convoluzione e sulla risposta impulsiva.

Un procedimento analogo può essere svolto nel caso continuo, de- Nel caso analogico èpossibile definire lederivate nonintere. . .

finendo la derivata come la convoluzione nel tempo tra la derivata diuna delta di Dirac ed il sistema [30, 28, 29]. Per la proprietà della deltadi Dirac [24, pagina 131] infatti

dδ(t)dt⊗ f (t) =

d f (t)dt

,

con la proprietà la proprietà della convoluzione (1.24) e della derivatadella trasformata di Laplace (1.26), si ottiene che

L

{dδ(t)

dt⊗ f (t)

}= L

{d f (t)

dt

}= sαF(s),


per cui antitrasformando ho la convoluzione tra la funzione f (t) e. . . comeconvoluzione per la

risposta impulsiva disα;

la risposta hα(t) del sistema sα, vale a dire la derivata di ordine αdell’impulso di Dirac

d f (t)dt

= hα(t)⊗ f (t),

che può essere considerata come une definizione di derivate fraziona-rie (in verità analoga alla definizione di Riemann-Liouville!).

Per sistemi discreti

Si può fare la stessa cosa nel caso discreto per ottenere una formulacon un procedimentoanalogo nel caso

discreto. . .valide nel caso discreto per le differenze frazionarie. Non si confonda ilsimbolo Dαx[n] con quello delle derivate Dαx(t) che si applica solo afunzioni continue.

Ricalcando il procedimento svolto nel caso analogico, per prima cosascriviamo le differenze frazionarie come convoluzione

Dαx[n] = hα[n]⊗ x[n], (3.11)

che possiamo riscrivere sostituendo il valore calcolato per la rispostaimpulsiva hα[n] (3.8) nella forma

Dαx[n] =sin[π(n− α)]

π(n− α)⊗ x[n],

che esplicita i termini della convoluzione in funzione dei campioniinteri del segnale.

Sappiamo in oltre che la differenza di ordine α è definita

aDαnx[n] = x[n− α],

traslando il segnale di una quantità α, per cui per confronto con la (3.8)possiamo assimilare, come a questo punto è ovvio che sia,. . . si ricava la

corrispondenza tra lepotenze e le differenze

di ordine α! z−αX(z)⇔ x[n− α], (3.12)

in modo del tutto analogo alla(3.5) [29]!In questo modo sono state generalizzate tutte le proprietà che ci

permettono di generalizzare l’equazione alle differenze finite di ordineintero

N

∑n=0

any[t− n] =M

∑m=0

bmx[t−m],

con l’equazione alle differenze non intere

N

∑n=0

any[t− αn] =M

∑m=0

bmx[t− αm], (3.2)

come avevamo già immaginato all’inizio di questa sezione.


3.2.2 Risposta in frequenza

È stato dimostrato nelle righe precedenti (3.12), che esiste una cor-rispondenza tra le potenze di z di ordine non intero α, e le diffe-renze finite traslate della stessa quantità α. Pertanto possiamo pas-sare dall’equazione alle differenze non intere (3.2), all’equazione infrequenza

N

∑n=0

ane−jωαnY(jω) =M

∑m=0

bme−jωαm X(jω),

dalla quale otteniamo infine

H(jω) = ∑Mm=0 bme−jωαm

∑Nn=0 ane−jωαn

, (3.13)

equazione che rappresenta la risposta in frequenza del circuito!

3.2.3 Trasformata Z

Finora siamo riusciti a ricavare le relazioni fondamentali per un siste-ma frazionario, solamente sulla circonferenza unitaria |z| = 1, per cuisiamo riusciti a ricavare la risposta in frequenza, ma non la funzione ditrasferimento del circuito.

Del resto, sono stati illustrati i motivi per i quali non è possibile L’antitrasformataZ −1 di z−α, insenso classico, non èdefinita. . .

sfruttare la definizione di antitrasformata Z ordinaria per calcolarela risposta di un termine della forma z−α, poiché tale termine nonè riconducibile ad alcuno sviluppo in serie di MacLaurin o Laurentche sia convergente, per cui non è possibile ricavare alcuna sequenzacausale o non causale, nel modo usuale.

Introduciamo quindi una definizione di trasfomata Z modificata,basata sulle formule integrali di Cauchy per proiettare una funzioneanalitica sulla circonferenza unitaria, su tutto il piano complesso, esten-dendo le relazioni valide solo per la DTFT.

Per le formule integrali di Cauchy [39], è possibile calcolare il valore . . . ma si può ottenereproiettando con gliintegrali di Cauchyla DTFT, definitasulla solacirconferenzaunitaria

assunto da una funzione f in un punto z come integrale sulla curva γche racchiude z

f−(z) =1

2π j

∮γ

f (w)w− z

dw, |z| < 1, (3.14a)

in cui f (z) è una funzione olomorfa4 in un dominio del piano complessoC in cui risiede la curva γ, percorsa in senso antiorario. Si è indicatocon f−(z) il sistema anti-causale, costituito cioè da sequenze unilateresinistre, la cui regione di convergenza è la proiezione all’interno delcerchio unitario che si ottiene con l’integrale di Cauchy.

Per sequenze unilatere-destre che danno origine a sistemi causali, l’in-tegrale di Cauchy assume la forma

f +(z) =1

2π j

∮γ

f (w)z− w

dw, |z| > 1, (3.14b)

in cui sono scambiati i termini al denominatore nell’argomento dell’in-tegrale.

4 Si ricorda che una funzione è olomorfa o analitica nel punto z = z0, se esiste un intorno|z− z0| < δ, in tutti i punti del quale esiste f ′(z) [39].


La funzione f (w) che siamo interessati a proiettare, è la rispostain frequenza di una sequenza x[n], che si calcola, come noto, dalladefinizione causale di Discrete Time Fourier Transform (DTFT):

X(ejω) ≡∞

∑k=0

x[k]e−jωk,

che tramite la sostituzione w→ ejω, mi permette di scrivere la funzionef (w), in forma più compatta

f (w) = X(w) =∞

∑k=0

x[k]wk. (3.15)

Gli integrali di Cauchy (3.14), calcolati sulla risposte in frequen-za (3.15), assumono quindi le forme

X−(z) =1

2π j

∮γ

∑∞k=0 x[k]wk

w− zdw, |z| < 1, (3.16a)

X+(z) =1

2π j

∮γ

∑∞k=0 x[k]wk

z− wdw, |z| > 1, (3.16b)

che possono essere assimilate alla definizione di trasformata Z modi-ficata! In particolare la definizione bilatera sarà data dalla somma delledue (3.16),

Xα(z) = X−(z) + X+(z),

in modo da risultare estesa all’intero piano complesso!

3.3 SISTEMI ELEMENTARIRicalcando il procedimento svolto al § 2.2 nella pagina 22 per i sistemitempo continui, cerchiamo di ricavare una procedura sistematica che cipermetta di affrontare e risolvere un sistema frazionario discreto [29].

Le basi poste fin qui, che qui riscriviamo, sono:Equazioni basilari deisistemi non interi

• l’equazione alle differenze finite non intere (3.2)

N

∑n=0

any[t− αn] =M

∑m=0

bmx[t− αm]; (3.2)

• la risposta in frequenza (3.13)

H(jω) = ∑Mm=0 bmz−jωαm

∑Nn=0 ane−jωαn

; (3.13)

• infine, le equazioni integrali di Cauchy (3.14) per estendere al pia-no complesso la risposta in frequenza definita sulla circonferenzaunitaria.

Oltre alle precedenti osservazioni, occorre prestare attenzione nel di-stinguere sistemi descritti da modelli Auto Regressive Moving Average(ARMA) frazionari (FARMA), cioè modelli a media mobile autoregressi-vi, o sistemi più semplici descritto solo, ad esempio, da un model-lo a media mobile. Queste distinzioni ci permettono di semplificarel’analisi [40, Cfr. Modelli stocastici].


3.3.1 Modelli MA

Nel caso di modelli a media mobile (Moving Average, ma), la risposta Per i modelli a mediamobile. . .in frequenza

H(ω) =M

∑m=0

bme−jωαm , (3.17)

è una forma polinomiale molto semplice, altrettanto semplice risulta ilcalcolo della risposta impulsiva h[n] e della funzione di trasferimentoH(z).

Sostituendo infatti la risposta in frequenza (3.17) nelle formule inte- . . . si semplifica ilcalcolo degli integralidi Cauchy

grali di Cauchy (3.14), otteniamo gli integrali

H+(z) =1

2π j

∫ π

−π

∑mm=0 bme−jωα

z− ejω jejωdω,

per il sistema causale, mentre

H−(z) =1

2π

∫ π

−π

∑mm=0 bme−jωα

1− ze−jω dω,

per il sistema anti-causale.Sviluppando le equazioni, ad esempio per H−(z), si ottiene la fun-

zione di trasferimento [29]

H−(z) =∞

∑n=0

z−nM

∑m=0

bm1

2π

∫ π

−πe−jω(n−m)αdω =,

e notando che l’integrale corrisponde alla derivata di ordine n − mαdell’impulso δ, cioè alla risposta impulsiva hn−mα di un sistema diordine frazionario n−mα (3.8), otteniamo

=∞

∑n=0

M

∑m=0

bmhn−mαz−n,

per cui, confrontando con la definizione di trasformata z, otteniamo la . . . si ottiene unarisposta impulsivapiuttosto semplice!

risposta impulsiva

h[n] =M

∑m=0

bmhn−mα. (3.18)

3.3.2 Modello completo

Nel caso di un modello completo Fractional Autoregressive MovingAverage (FARMA), la risposta in frequenza completa da considerare èla (3.13), con N > M per avere un sistema proprio.

Anche nel caso discreto è possibile ricondurre gli esponenti αi5 a Esprimo gli

esponenti comemultipli interi di unsolo esponentefrazionario.

multipli di un valore base, per cui nel considerare la(3.13), è possibiledistinguere i casi:

1. Gli esponenti αi sono numeri razionali. In tal caso, dopo aver tro-vato il minimo comune multiplo α, indichiamo tutti gli esponenticome suoi multipli, iα per i = 0, 1, 2, . . . , I;

5 in questo caso con i si intendono sia gli indici n che m, con I invece si intende sia M cheN.


2. Gli esponenti α sono numeri reali, ma multipli di un valore dibase α ∈ (0, 1). Ancora una volta è possibile scrivere gli espo-nenti come multipli di α, quindi si ottiene comunque αi = iα peri = 0, 1, 2, . . . , I.

In definitiva in entrambi i casi la trattazione è la stessa, si parte dauna funzione di trasferimento che è stata posta nella forma

H(jω) = ∑Mm=0 bme−jωmα

∑Nn=0 ane−jωnα

,

poi si seguono i seguenti semplici passi:

1. trasformiamo H(z) in H(w), sostituendo zα con w. Se H(w)non fosse una frazione propria, sarebbe necessario estrarre laparte impropria come somma polinomiale, analogamente al casocontinuo;

2. effettuare un’espansione in frazioni parziali della H(z), per otte-nere una funzione data dalla somma di termini elementari;

3. torniamo nuovamente in z, con la sostituzione di zα al posto w;

4. invertiamo ciascuna funzione elementare e sommiamo le risposteimpulsive ottenute, per ottenere la risposta impulsiva complessi-va del sistema.

Le proprietà della risposta impulsiva, sono ereditate dalle proprietàIl sistema è stabile sei poli sono esterni al

cerchio unitariocorrispondenti a ciascun fratto semplice, per cui è possibile studiare lastabilità del sistema dalla posizione dei poli di ciascun fratto semplice.Come nel caso intero, se i poli sono tutti all’esterno della circonferenzaunitaria, il circuito risulterà stabile.

3.3.3 Inversione fratti semplici

Resta come ultimo problema, capire come sia possibile invertire cia-scuno dei termini in cui abbiamo scomposto la funzione di trasferi-mento [29]. Consideriamo quindi un termine semplice del tipo:

F(w) =1

1− aw−α(3.19)

definito per w = e−jω.Riprendendo gli integrali di Cauchy (3.14) per l’inversione ottenia-

mo

F(z) =1

2π j

∮γ

11− aw−α

1z− w

dw, |z| > 1,

per i sistemi causali, da cui

F(z) =1

2π j

∮γ

11− ae−jωα

1z− ejω dw, |z| > 1,

sfruttando la condizione w = e−jω che il dominio sia sulla circonferen-za unitaria!


Una soluzione semplice del problema si ottiene elaborando l’inte-grale in modo che assuma la forma di una trasformata. In pratica perprima cosa sfruttiamo l’equivalenza

11− ae−jωα

=∑

q−1k=0 ejkωα

1− aqe−jω , (3.20)

valida nel caso in cui α = 1/q, come si dimostra facilmente calcolandoil valore assunto dalla sommatoria al 2° membro.

Inoltre, se |a| < 1, si può riscrivere la (3.20) nella forma

11− ae−jωα

=∞

∑m=0

q−1

∑k=0

aqm+ke−jω(m+kα),

per cui possiamo infine calcolare la trasformata Z inversa di uno deifratti semplici (3.19) in cui abbiamo scomposto la funzione, ottenendoinfine la risposta impulsiva della forma:

h[n] =q−1

∑k=0

ak∞

∑m=0

aqmhnm−kα (3.21)

notando anche in questo caso, che l’integrale corrisponde alla derivatadi ordine nm− kα dell’impulso δ, cioè alla risposta impulsiva hnm−kα

di un sistema di ordine frazionario nm− kα (3.8).Come ultima considerazione, osserviamo infine che mentre per i mo-

delli a media mobile la risposta impulsiva del sistema è finita (3.18),per i modelli FARMA completi otteniamo una risposta impulsiva (3.21)che necessita di uno sviluppo illimitato, saranno necessarie quindidelle approssimazioni per poter avere un modello implementabile.

3.4 APPROCCIO ALTERNATIVOPer quanto riguarda i sistemi discreti di ordine frazionario, come mo-strato nelle sezioni precedenti, esistono numerosi approcci possibili.Tutti i metodi si basano sull’approssimazione della funzione di trasfe-rimento del sistema; come funzione di trasferimento di ordine intero Sistemi discreti non

interi: assenza dimetodi alle differenzefinite

(con un numero più elevato di poli e zeri), o come funzione di z cheinterpoli il segnale tra i suoi campioni. Non è stato possibile trovarein nessun articolo un’approssimazione alle differenze finite, che facciauso dei soli campioni del segnale.

Cerchiamo quindi di trovare una relazione tra i campioni dei segnalidi ingresso e di uscita, che si possa interpretare come una derivatadi ordine non intero nel dominio continuo, e che possibilmente siasufficientemente pratica ed accurata nei risultati.

Osserviamo 2 equazioni ben note nell’analisi di sistemi interi, la pri- Analogia con isistemi discreti noninteri

ma definita in un dominio continuo, l’altra, equivalente alla prima, neldominio discreto:

N

∑k=0

akDky(t) =M

∑k=0

bkDkx(t), (2.1)

N

∑k=0

ak y[n− k] =M

∑k=0

bk x[n− k], (3.1)


in cui le derivate della prima sono approssimate da differenze finitenella seconda. Si vuole ricavare un’equazione analoga all’equazionealle differenze (3.1), capace di approssimare un’equazione differenzialedi ordine non intero!

3.4.1 Discretizzazione delle derivate non intere

Consideriamo quindi un sistema tempo-continuo non intero; l’equazio-ne (2.1) scritta per k interi, si scrive, come mostrato nell’equazione (2.2),come la sommatoria di N termini derivati che indicheremo con αk.Esplicitiamo la derivata con la definizione di Grünwald-Letnikov (1.13)Dalla definizione di

G-L data nel § 1.2.2 nella pagina 4, che qui riportiamo

aDαt x(t) =

(t− a

N

)α N−1

∑k=0

Γ(k− α)Γ(−α)Γ(k + 1)

x(

t− kt− a

N

)(3.22)

con una prima modifica: è stato rimosso il limite per N → ∞ presentenella definizione originale. Questa prima approssimazione ci permettedi trattare con segnali x ed y campionati, cioè noti solo negli istantidi tempo t = nT; il numero di campioni che è possibile gestire inqualsiasi applicazione reale deve essere necessariamente finito.

Osserviamo poi, che il differenziale quantizzato, che nella formu-Il differenzialequantizzato come

periodola di G-L è scritto come (t − a)/N, equivale al periodo T del segnalecampionato N + 1 volte tra un estremo inferiore fisso a e l’estremosuperiore t. Lo stesso concetto in formule:

t− aN

= T,

sostituito nella (3.22) ci permette di scrivere:

aDαt x(t) = Tα

N−1

∑k=0

Γ(k− α)Γ(−α)Γ(k + 1)

x(t− kT)

Il rapporto tra le tre funzioni Γ non dipende dall’istante t, ma soloDipendenza dellasola x dal tempo t dall’indice della sommatoria k corrente e dall’ordine di derivazione α.

È possibile pensare quindi di sfruttare una notazione più compatta,in cui i tre termini ed il fattore pari ad una potenza di α del periodo,siano inglobati in un coefficiente, che posso indicare banalmente conΛα[k], ottenendo

aDαt x(t) =

N−1

∑k=0

Λα(k) x(t− kT)

Considerando infine di operare su segnali campionati, è possibile rife-Equivalenza tratempo e campioni rirsi direttamente all’indice numerico piuttosto che all’istante di tempo

in cui campiono: è possibile quindi scrivere la sostituzione: t− kT →n− k, ottenendo infine

Dαx[n] =N−1

∑k=0

Λα[k] x[n− k], (3.23)

in cui le parentesi quadre indicano che l’argomento della funzione èdiventato un indice del campione considerato.

Un modo ancora più compatto per indicare l’n-esimo campione del-Notazione vettoriale


la derivata di ordine α della funzione si ottiene con la notazione vetto-riale. L’argomento della funzione x, per k = 0, coincide con l’estremodestro t dell’intervallo considerato; per k crescenti, decresce di k perio-di T, cioè di k campioni del segnale. Indicando sia il segnale x, sia icoefficienti Λα[k], con un vettore colonna con le righe ordinate per kcrescenti, cioè:

Λα ≡ Tα ·

Λ[0]Λ[1]Λ[2]

...Λ[N − 1]

, x ≡

x[N − 1]x[N − 2]x[N − 3]

...x[0]

,

si arriva alla notazione:

Dαx[n] = ΛTα · x, (3.24)

estremamente compatta che sfrutta il prodotto righe per colonne dellematrici. Ovviamente la T indica che il vettore Λα deve essere trasposto.

3.4.2 Sistemi descritti da equazioni alle differenze finite

Torniamo adesso alla corrispondenza tra i sistemi tempo continui etempo discreti, che per ordini interi sono descritti rispettivamente dalleequazioni (2.1) ed (3.1). La stessa analogia per sistemi non interi si Equazione vettoriale

alle differenze finiteottiene tramite la (3.24) a partire dall’equazione differenziale (2.2) perun sistema tempo continuo di ordine non intero. Si ottiene facilmentela seguente

N

∑n=0

an ΛTαn · y =

M

∑m=0

bm ΛTαm · x (3.25)

relazione vettoriale tra i segnali di ingresso e di uscita.Per esplicitare il ruolo dei ritardi, mascherati dalla notazione vetto- Equazione esplicita

alle differenze finiteriale, si può sostituire l’equazione (3.23) nella (3.1), ottenendo

N

∑n=0

an

C−1

∑k=0

Λαn [k] y[n− k] =M

∑m=0

bm

C−1

∑k=0

Λαn [k] x[n− k], (3.26)

che qui confrontiamo riportando la corrispondente equazione conti-nua(2.1):

N

∑n=0

anDαn y(t) =M

∑m=0

bmDαm x(t).

L’equazione (3.26) permette di esprimere la relazione tra l’ingresso e Considerazionil’uscita di un sistema in termini degli ingressi e delle uscite presenti epassate, analogamente ad una comune equazione alle differenze finitedi un sistema intero. Non è stato possibile trovare in letteratura descri-zioni di sistemi che siano assimilabili all’equazione data. Molti degliaspetti relativi all’approssimazione delle derivate frazionarie, applica-ti per ricavare tale equazione, si trovano sui libri di analisi non intera,come ad esempio [27].


3.4.3 Algoritmi per equazioni alle differenze

L’equazione (3.26) è analoga per molti aspetti ad un’equazione persistemi a coefficienti interi come la (3.1).

Una prima differenza è la presenza di un numero elevato di ter-Sistemi discreti noninteri operano con un

numero infinito ditermini troncato

mini passati, per C punti impiegati nel calcolo ho C termini passati,anche nel caso di una sola derivata. Per sistemi interi invece, se hosolo una derivata, ad esempio per M = 1 nella (2.1), ho un sistemadiscreto la cui uscita y che dipende solo dall’ingresso corrente x[n] edal precedente x[n− 1].

Del resto, riducendo il numero di termini C coinvolti nel calcolo,ho un’approssimazione peggiore. Proviamo a valutare quanto inci-de il numero di punti C sull’accuratezza. Consideriamo la funzioneMatlab:

1 function G=vg(q,N)

che prende in ingresso l’ordine di derivazione α indicato nel codicecon q, ed il numero di punti C, nel codice indicati con N, per restituireil vettore:

Λα ≡

Λ0Λ1Λ2...

ΛN−1

, Λα[k] =Γ(k− α)

Γ(−α)Γ(k + 1)per k = 0, 1, . . . , N− 1,

(3.27)

vettore delle tre funzioni gamma, i cui termini dipendono esclusiva-mente dall’ordine di derivazione. La dipendenza da k è fittizia, nelsenso che limitare i k possibili ai primi N − 1 valori interi più lo zero,indica l’approssimazione applicata, non un reale grado di libertà. Ilvettore inoltre è stato normalizzato rispetto al termine costante T−α.L’algoritmo che calcola gli elementi del vettore di equazioni GammaCalcolo del vettore

Gamma definito dalla (3.27), è dato dalle semplici istruzioni seguenti.

3 G=1;

4 for i=1:N-1

5 G(i+1)= G(i) * (i-1-q)/i;

6 end

Listato 3.1: Algoritmo di calcolo del valori dei coefficienti funzioni Gamma

Analisi del codice

Il codice sfrutta la proprietà ricorsiva della funzione Gamma. Esplici-tando i primi elementi è possibile scrivere:

Λ =1

Γ(−α)·{

Γ(−α)Γ(1)

,Γ(1− α)

Γ(2),

Γ(2− α)Γ(3)

, . . .}

,

ma per la proprietà ricorsiva (1.4), il 2° termine si può esprimere infunzione del 1°,

Γ(1− α)Γ(2)

=−αΓ(−α)

1Γ(1)= −α Λα[0],


il 3° termine in funzione del 2°, quindi del 1°,

Γ(2− α)Γ(3)

=(1− α)Γ(1− α)

2Γ(2)=

1− α

2·Λα[1] =

−α(1− α)2

·Λα[0],

e così via. In generale quindi, il termine Γ(−α) si semplifica, la ricor- formula ricorsiva peril calcolo dellafunzione Gamma

sione si può scrivere:

Γ(k− α)Γ(k + 1)

=k− 1− α

kΓ(k− 1− α)

Γ(k), (3.28)

abbiamo così ottenuto la funzione implementata nel codice 3.1.

Considerazioni

Il metodo proposto in questa sezione, è stato pensato con l’intento di Metodi ‘nativi’ alledifferenze non interedefinire un’equazione per i sistemi non interi analoga alla (3.1). In let-

teratura sono disponibili, come riportato in questa sezione, approfon-dimenti sul metodo di studio di sistemi frazionari, concepiti già comesistemi discreti, ovvero concepiti come sistemi di equazioni alle diffe-renze non intere. Questi non si ottengono come discretizzazione disistemi analogici, ma possono essere concepiti nella forma (3.2) a parti-re da considerazioni matematiche sul sistema; un esempio applicativosarà considerato nel § 4.4.2 nella pagina 80.

L’equazione (3.25) propone invece un’espressione compatta per i si-stemi discreti di ordine non intero, che deriva direttamente da un’e-quazione differintegrale. Al tempo della stesura, non è stata trovataalcuna fonte che proponesse qualcosa di simile. Successivamente l’ap-proccio vettoriale proposto è stato ritrovato in [16], che sviluppa la In [16] è proposto un

metodo matematicodi calcolo cheequivale al metodoqui presentato

parte matematica in modo analogo, sebbene tratti l’argomento esclusi-vamente come metodo di calcolo, tralasciando completamente tutte lealtre tematiche relative ai sistemi discreti.

Matlab dispone di una funzione per il calcolo della funzione gam- La funzione gammain Matlabma, che prende come argomento un qualsiasi numero reale. Ricor-

dando che la funzione gamma cresce come il fattoriale per ingressipositivi, ci si imbatte presto nella limitatezza del calcolatore usato, in-fatti l’istruzione gamma(171) restituisce 7.2574e+306, mentre l’istruzionegamma(172) restituisce il valore Inf, per cui non è possibile calcolare lafunzione gamma per ingressi maggiori di 172!

Con l’algoritmo mostrato nella (3.28) non è necessario calcolare espli-citamente la funzione gamma, per cui posso calcolare vettori anche conpiù di 171 elementi senza incappare nei limiti considerati.

D’altro canto, si osserva inoltre che un numero di elementi elevatospesso non è necessario in pratica, poiché il vettore Λα decresce con-siderevolmente per k crescenti, come mostrato nel grafico 5: nel casopeggiore bastano una dozzina di elementi affinché il peso attribuitoalla funzione scenda al di sotto di un centesimo del peso del primoelemento. Necessità di una

formula ricorsivaCome vedremo nel § 4.3 nella pagina 70, la formula ricorsiva quitrovata per il calcolo della funzione gamma, può essere estesa al cal-colo generale di un coefficiente binomiale con termini non interi inmodo efficienti. Il procedimento richiede che i coefficienti binomialisiano riscritti con la (1.6), valida anche per n negativi; successivamentesi sostituisce al fattoriale la funzione gamma (1.3) ottenendo una for-mula analoga alla (1.7); infine si riscrive la formula ottenuta in modoricorsivo tramite la (1.4).


0 5 10 15 20 25−35

−30

−25

−20

−15

−10

−5

0

# Elemento del vettore

Valo

re (d

B)

Andamento del vettore Gamma

k=0.50.850.10.251

Figura 5: Andamento degli elementi del vettore Λα[k] al variare del parametroα. Si osservi che l’asse delle ordinate è in dB!

I tentativi svolti in Matlab, per il calcolo di un coefficiente bino-miale tramite la definizione diretta, non hanno mai portato a risulta-ti ammissibili, anche per coefficienti binomiali piccoli, come (0.3

5 ), cherientrano comodamente nei limiti consentiti dal calcolatore in uso.

3.5 APPROSSIMAZIONI DISCRETEAnalogamente a quanto fatto nel dominio tempo continuo, esistono va-rie tecniche per approssimare un sistema discreto. In generali i metodipossibili si possono raccogliere in due gruppi:

• Discretizzazione diretta;

• Discretizzazione indiretta.

Nei metodi indiretti, si hanno due passaggi: per prima cosa si trasfor-Metodi indirettima il filtro tempo continuo in frequenza e si approssima per avere unfiltro nella variabile di Laplace s di ordine intero, successivamente sidiscretizza la variabile di Laplace con una trasformazione. In praticaoccorre applicare uno dei metodi descritti nel § 2.4 nella pagina 27, epoi applicare un metodo classico di discretizzazione.

Nei metodi diretti invece, dei quali ci occuperemo in questa sezione,Metodi Direttisi discretizza direttamente la variabile di Laplace di ordine non inte-ro, ad esempio tramite CFE della trasformazione di Eulero, o tramitemetodi di integrazione numerica.

3.5.1 Necessità delle approssimazioni

Un sistema digitale richiede una equazione alle differenze finite, chenel caso intero corrisponde ad una equazione differenziale nel domi-nio continuo. Per i sistemi frazionari non si ha una corrispondenzaaltrettanto chiara tra i domini continuo e discreto [30]; campionandocon un periodo T si ottiene una equazione alle differenze finite basata


su ritardi frazionari, che quindi può essere messa in relazione con unsistema ordinario con decimazione o interpolazione di campioni!

Sono stati studiati più metodi di approssimazione di equazioni, perpassare dal caso continuo a quello discreto nei domini trasformati, cioèdalla trasformata di Laplace alla trasformata Z , e per approssima-re funzioni di trasferimento frazionarie con funzioni di ordine intero.Vediamo nelle prossime sezioni quali sono le possibili operazioni diapprossimazione e le relative peculiarità.

Le funzioni di discretizzazione che analizzeremo permettono di pas- discretizzazione dallarisposta nel dominiodi Laplace

sare dalla variabile di Laplace sα ad un’equazione in z equivalente chel’approssima secondo determinati criteri. Le approssimazioni discre-te di Eulero Grünwald-Letnikov , di Al-Alaoui e l’approssimazionebilineare di Tustin stabiliscono una corrispondenza tra le trasforma-te L e Z , cioé tra la variabile di Laplace sα ed una funzione di z,analogamente a quanto avviene per funzioni di ordine intero.

Ripartiamo dalle equazioni continue del sistema, nel dominio deltempo come la (2.2), o nel rapporto ingresso/uscita nel dominio diLaplace, cioè nella forma di funzione di trasferimento (2.4). Per ot-tenere un modello discreto è necessario approssimare derivate ed in-tegrali frazionari, cioè i termini sα, con un loro equivalente discreto.L’approssimazione che si cerca dovrebbe essere: Caratteristiche

dell’approssimazionecercata• in forma frazionale, con poli e zeri, che dovranno essere inter-

lacciati per approssimare pendenze generiche con multipli di20 dB/decade;

• sufficientemente accurata anche con un ridotto set di parametri,o sufficientemente facile da calcolare per poter essere implemen-tabile in tempo reale;

• il semipiano sinistro della variabile complessa di Laplace dovreb-be essere mappato nella circonferenza unitaria del piano Z , conl’asse immaginario che corrisponde alla circonferenza stessa [2].

In generale, si può esprimere l’approssimazione z con una funzio-ne w(z), detta funzione generatrice. La (2.4), diventa quindi, persostituzione della variabile di Laplace s con un’opportuna funzionegeneratrice w(z),

H(z) = ∑Mm=0 bm[w(z)]αm

∑Nn=0 an[w(z)]αn

.

Le funzioni generatrici tipicamente utilizzate sono le stesse applica- Metodi delle funzionigeneratricibili al caso intero, elevate però ad esponente non intero, cioè

• metodo di Eulero Grünvald-Letnikov;

• metodo di Tustin o trapezoidale (bilineare);

• metodo di Al-Alaoui.

ciascuna delle quali corrisponde ad un’espressione frazionale in z.Ciascuna approssimazione è illimitata, perciò richiede una ulteriore Funzioni non

limitate erealizzabilità

approssimazione, dalla quale è possibile ricavare la serie dei coeffi-cienti della risposta impulsiva, FIR nel caso di espansione polinomialetroncata, IIR nel caso di espansione CFE. Infatti, se consideriamo il me-todo di Eulero, tanto per fissare le idee, questo si basa sulle differen-ze finite all’indietro 1− z−1, infatti dalla sua espansione Power Series


Expansion (PSE) si ottiene la definizione di differintegrale di Grünwald-Letnikov. Elevando a potenza intera n ottengo un polinomio di gradon, ma se l’esponente è frazionario, ottengo una serie infinita di termini,proprio come avviene per le derivate non intere (1.13) di G-L.

Approssimazioni della (2.2), si possono invece ottenere tramite al-goritmi di integrazione numerica, in letteratura tali metodi sono pocodiffusi ed illustrati, per questo si rimanda alla sezione 3.4 nella pa-gina 49. Lo svantaggio di tali metodi è nelle maggiori richieste dimemoria.

Analizziamo ora nelle prossime pagine i vari metodi di approssi-mazione disponibili, che oltre a quelli citati comprendono i metodiderivati dalle tecniche di elaborazione dei segnali di Padé, Prony eShanks, dai quali si ottengono filtri IIR.

3.5.2 Metodo dell’errore quadratico minimo

Abbiamo visto, (3.2) nella pagina 40, che un sistema discreto non in-tero presenta un’equazione alle differenze non intere, e che per cia-scun termine della (3.2) si ha una risposta impulsiva di tipo sinc, comemostrato nella (3.10) nella pagina 43.

Il ricampionamento ideale, con la risposta impulsiva di tipo sincrichiederebbe infiniti termini, come già detto, ciò non è possibile inpratica, per cui sfruttiamo un’approssimazione finita di L termini, chesi avvicini il più possibile alla risposta ideale.

Per fissare le idee, supponiamo di avere un segnale x[n], elaboratoin un sistema non intero con ritardo α, in un’equazione ingresso uscitatipo y[n] = x[n − α]. Consideriamo il segnale x[n − α] traslato di α,come l’uscita del segnale x[n] da un sistema che introduce un ritardoα. Come visto al §3.2, possiamo scrivere

xi[n− α] = x[n]⊗ hα[n] =+∞

∑k=−∞

x[k]hα[n− k],

in cui hα[n] è proprio la risposta impulsiva del sistema ritardante otte-nuta nella (3.10), e la i a pedice della x ci ricorda che stiamo conside-rando un caso ideale.

Passando al dominio Z , si può scrivere [23, 34]

Hi(z) =+∞

∑n=−∞

hα[n]z−n

per la funzione di trasferimento ideale, e

H(z) =b

∑n=a

hα[n]z−n

per la funzione di trasferimento approssimata con L termini (L = b−a + 1).

Il minimo dell’errore quadratico6, si ottiene minimizzando la norma 2

dell’errore

E(jω) = H(ejω)− Hi(ejω)

6 I metodi dell’errore quadratico sono spesso conosciuti con il nome di least square, termineusato di frequente anche nei testi italiani.


in un intervallo di frequenze [0, π], per cui si avrà

‖E‖2 =1pi

∫ π

0|E(ejω)|2dω =

1π

∫ π

0|H(ejω)− Hi(ejω)|2dω. (3.29)

Grazie al teorema di Parseval che ricordiamo:∫ +∞

−∞|x(t)|2dt =

∫ +∞

−∞|X( f )|2d f ,

possiamo esprimere la (3.29) in funzione delle risposte impulsive neltempo

e[n] =1π

−∞

∑−∞|h[n]− hi[n]|2,

dalla quale possiamo dedurre che otteniamo il minimo dell’errore sem-plicemente troncando la riposta impulsiva, considerando solo gli L ter-mini centrati attorno all’origine della funzione sinc, e che, ovviamente,l’accuratezza è tanto maggiore quanti più termini consideriamo.

Gli estremi a, b, sono quindi presi alla stessa distanza dal centro,cioè

a =

{[α]− L

2 , per L pari,bαc − L−1

2 , per L dispari,

in cui bαc è la parte intera di α, cioè il suo valore approssimato perdifetto, mentre [α] è l’approssimazione intera più prossima al suo va-lore, per eccesso o per difetto. Il filtro quindi, sarà nullo all’esternodell’intervallo [a, a + L], e l’errore residuo risulterà

‖E‖2 =a−1

∑n=−∞

|hi|2 +∞

∑n=a+L+1

|hi|2,

mentre l’equazione di un sistema della forma (3.2), con l’approssima-zione proposta si esprime

N

∑n=0

an

L

∑k=a

y[k] sinc[π(n− k− αn)] =M

∑m=0

bm

L

∑k=a

x[k] sinc[π(n− k− αm)],

con L termini di approssimazione.

3.5.3 Metodo di Eulero-Grünwald

Nel metodo di Eulero frazionale si approssima l’equazione differenzialedel sistema (2.1), con le differenze finite all’indietro7, passando poi aldominio Z , per ottenere la trasformazione Funzione generatrice

del metodo

s ≈ 1− z−1

T, (3.30)

ben nota per i sistemi di ordine intero, anche detta metodo delle diffe-renze finite all’indietro.

7 Ovviamente il caso di differenze in avanti, più rigoroso dal punto di vista matematico,non è applicabile per sistemi causali.


−4 −3 −2 −1 0 1 2 3 4−2

−1

0

1

2

3

4

5

x

1/Γ(x)

Figura 6: Funzione gamma inversa 1/Γ(x). Mentre per valori positivi la fun-zione gamma cresce come il fattoriale, quindi la sua inversa tendea 0, per valori interi negativi la funzione gamma tende ad infinito, ela sua inversa è nulla.

L’estensione al caso non intero è immediata, basta infatti sostituirela (3.30) nei sistemi del tipo (2.2) di ordine non intero, per ottenere [7]:

sα ≈[

1T

(1− z−1

)]α

. (3.31)

Cerchiamo di sviluppare il termine (1− z−1)α, tenendo in mente cheuna qualsiasi potenza intera n di un binomio, si può scrivere, come èfacile ricavare, con l’espressione

(a + b)n =n

∑k=0

(nk

)akbn−k, n ∈N

che si estende per n negativi, ponendo (nk) = 0 sia per n < 0 che per

n < k

(a + b)n =∞

∑k=0

(nk

)akbn−k, n ∈ Z

come mostrato, ad esempio, da [45].A questo punto possiamo sviluppare l’equazione (3.31), esplicitando

le potenze di binomio di ordine α tramite la proprietà fattoriale (1.3)della funzione gamma, ottenendo

(1− z−1

)α=

∞

∑k=0

(α

k

)(−z−1

)k=

∞

∑k=0

(−1)k(

α

k

)z−k, (3.32)

cioè una serie di potenze di ordine intero illimitata. Infatti, nel casonon intero non si verifica più la condizione (α

k) = 0 per α < k, come sipuò osservare dal reciproco della funzione gamma 6.


Osserviamo che i termini della serie infinita che si ottengono per α diordine non intero, coincidono con gli stessi delle derivate di Grünwald-Letnikov (1.13) o (1.11). In letteratura l’approssimazione

sα ≈(

1T

)α ∞

∑k=0

(−1)k(

α

k

)z−k

risultante, non è riferita allo sviluppo dei coefficienti binomiali (3.32),ma è riferita come Power Series Expansion (PSE), espansione in serie dipotenze, della (3.31).

La precedente equazione, raccogliendo i termini nella sommatoria,assume la forma di trasformata Z ,

sα ≈∞

∑k=0

(1T

)α

(−1)k(

α

k

)z−k

in cui emergono evidenti i coefficienti della risposta impulsiva hα[k] Risposta impulsivadel sistema sα campionato:

hα[n] =(

1T

)α

(−1)n(

α

n

)come mostrato da [7].

In tale risposta impulsiva non si ha alcuna retroazione dall’uscita,per cui è possibile ottenere un’implementazione FIR troncando lo svi-luppo ottenuto! Si ricorda però, che ai sistemi FIR nel caso non intero,corrisponde una scarsa accuratezza del modello analogico approssima-to.

Con un’espansione CFE invece, si può ottenere una funzione di tra-sferimento frazionale che si implementa con un filtro IIR [5], che ingenerale sarà più accurata

Con il metodo di Eulero, il semipiano s sinistro si mappa nella cir-conferenza con centro in z = 1/2 e raggio 1/2, per cui l’approssima-zione è abbastanza accurata solo in bassa frequenza.

Confronto con il metodo di differintegrazione numerica

Si è osservato, che la soluzione qui mostrata è analoga alla soluzionemostrata nel § 3.4 nella pagina 49. Infatti gli elementi della rispostaimpulsiva hα[n] possono essere espressi, grazie alla (1.7), nel rapporto(

1T

)α

(−1)n(

α

n

)= T−α Γ(k− α)

Γ(−α)Γ(k + 1)= Λα[k],

e quindi sostanzialmente i due metodi coincidono.Nella tabella 2 nella pagina successiva sono elencati i passaggi chia-

ve di entrambi i procedimenti, il risultato finale è lo stesso sebbenei percorsi siano concettualmente differenti. Nel caso della soluzionenumerica dei differintegrali, il vettore Λ coincide con la risposta im-pulsiva hα[n] del filtro frazionario del metodo di Eulero. Il prodottovettoriale esplicita la convoluzione in modo del tutto equivalente!

3.5.4 Metodo di Tustin

L’approssimazione di Tustin, o bilineare, o metodo trapezoidale, si ot-tiene trovando la soluzione del sistema tempo-continuo, approssiman-


Tabella 2: Confronto tra il metodo dell’approssimazione numerica deidifferintegrali di G-L, ed il metodo di Eulero frazionario.

Differintegrali di G-L Eulero Frazionario

Equazione integrodifferenziale nel tempoN

∑k=0

akDky(t) =M

∑k=0

bkDkx(t) (2.1)

Approssimazione delle derivatefrazionarie con la definizione diG-L, per un numero finito dicampioni della funzione:

Trasformazione nel dominio diLaplace della funzione con laregola di derivazione (1.26):

N

∑n=0

an ΛTαn · y =

M

∑m=0

bm ΛTαm · x

(3.25)

N

∑n=0

ansαn Y =M

∑m=0

bmsαm X (2.3)

Discretizzazione di s (3.31) esviluppo binomiale di ordineα (3.32) (in [n]):

N

∑n=0

anhα ⊗ y =M

∑m=0

bmhα ⊗ x

do l’integrale con il metodo trapezoidale. Come noto, la trasformazio-ne caratteristica del metodoFunzione generatrice

del metodo

s =2T· 1− z−1

1 + z−1 , (3.33)

è una trasformazione bilineare, che mappa il semipiano sinistro di snella parte interna del cerchio unitario del piano z.

Nel caso non intero, l’estensione [13]

s =(

2T· 1− z−1

1 + z−1

)α

, (3.34)

richiede solo l’aggiunta di un esponente frazionario.I coefficienti della risposta impulsiva si ottengono anche in questo

caso con uno sviluppo PSE, come mostrato da [7], si ottiene la rispostaimpulsiva

hα[n] =(

2T

)α n

∑k=0

(−1)k(

α

n

)(−α

n− k

), (3.35)

i cui coefficienti si possono calcolare sostituendo la funzione gammanello sviluppo fattoriale di coefficienti binomiali (1.7). La risposta im-pulsiva è illimitata come sempre avviene per esponenti non interi, perquesto deve essere troncata. I coefficienti possono essere efficientemen-te calcolati con un algoritmo ricorsivo per la funzione gamma, sfruttan-do la proprietà (1.4) fino a ricondursi alla funzione Γ(1), analogamentea quanto fatto al § 3.4.3 nella pagina 52.

Algoritmo ricorsivo di Chen e Moore [13]

L’approssimazione FIR della risposta impulsiva, è facilmente superatada un’approssimazione IIR, ottenuta ad esempio da uno sviluppo CFE


1 function B=Tustin(a,L)

2 if nargin<1,

3 L=9;

4 end

5

6 % Simbolico

7 z=sym(’z’);

8

9 A=1-a*z^(-1); % A1(z^-1, a)

10 B=1-a*z; % A1(z,a)

11 for i=2:L, % sono partito da A0,A1,continuo da A2 fino ad AL

12 if rem(i,2) %rem(i,2) = 1 per i dispari, quindi eseguo l’if

13 c=a/i;

14 else

15 c=0;

16 end

17 A = A - (c * B) * z^(-i);

18 B = subs(A,z^(-1));

19 beep

20 endListato 3.2: Codice per l’implementazione simbolica dell’algoritmo di Chen eMoore [13].

della (3.34). Un metodo conveniente per ricavare un’approssimazionefrazionale della (3.34) è stato proposto da Chen e Moore [13], basato Algoritmo ricorsivo

di Chen e Moore [13]su un algoritmo ricorsivo, per esprimere la funzione generatrice comerapporto di polinomi interi

(w(z)

)α =(

2T· 1− z−1

1 + z−1

)α

=(

2T

)α

limn→∞

An(α, z−1)An(−α, z−1)

. (3.36)

Il metodo ricorsivo permette di calcolare facilmente i polinomi An,sfruttando le relazioni di ricorrenza

An(α, z−1) = An−1(α, z−1)− cnz−n An−1(α, z),

cn =

{α/n per n dispari,0 per n pari,

a partire dalla condizione iniziale

A0(α, z) = 1.

Si osservi la presenza di due diversi polinomi A nell’algoritmo, che sidifferenziano per il segno dell’esponente di z, che si possono ottenereper semplice sostituzione delle variabili z↔ z−1.

Implementazione

L’algoritmo di Chen e Moore [13] è stato implementato in Matlab, conuna funzione che restituisce il valore simbolico, cioè un polinomio nellavariabile z, dell’approssimazione di ordine n del polinomio An.

Per semplicità, il polinomio restituito dallo script 3.2 è nella variabilez, i coefficienti sono gli stessi del polinomio nella variabile z−1.


3.5.5 Metodo di Al-Alaoui

Il metodo di Al-Alaoui [1, 2, 13] consiste in un’interpolazione dei me-todi di Eulero e Tustin, secondo un peso che può essere fissato a priorioppure variabile. Seguendo il metodo originale, l’interpolazione si ot-La funzione

generatrice è lasomma pesata dei

metodi di Eulero eTustin

tiene con pesi pari a 3/4 per la funzione di Eulero (HE) e 1/4 per lafunzione di Tustin (HT)

HA =34

HE +14

HT , (3.37)

cioè i pesi sono l’uno il triplo dell’altro, come mostrato nella formulaprecedente.

Sviluppando la (3.37) con le espressioni esplicite delle funzioni diEulero (3.30) e di Tustin (3.33), otteniamo la

sα ≈(

87T

1− z−1

1 + z−1

7

)α

,

formula generatrice esplicita del metodo di Al-Alaoui.Alcuni articoli ([15, 14, 25]) suggeriscono l’uso della funzione gene-

ratrice di Simpson

s =[

3T· (1 + z−1)(1− z−1)

1 + 4z−1 + z−2

]α

con un peso variabile, al posto della funzione di Eulero, definendocomunque come di Al-Alaoui la funzione risultante.

Anche in questo caso è possibile ottenere i coefficienti della rispostaimpulsiva da uno sviluppo PSE

hα[n] =(

87T

)α n

∑k=0

(−1)k(

17

)n−k(α

n

)(−α

n− k

)(3.38)

La funzione di Al-Alaoui è la più accurata nel tracciamento di fun-zioni in alta frequenza, anche in questo caso è possibile un’espansioneCFE per ottenere una forma FIR realizzabile!

3.5.6 Padé, Prony e Shanks

Sono stati sviluppati degli efficaci metodi di approssimazione derivatidalle tecniche di elaborazione dei segnali, proposti da Padé, Prony eShanks [7, 6], che permettono di avere un’approssimazione IIR, cioèdella forma

H(z) =B(z)A(z)

=b0 + b1z−1 + · · ·+ bmz−m

a0 + a1z−1 + · · ·+ anz−n , (3.39)

con m ≤ n affinché sia una frazione propria.Il filtro viene determinato quando vengono fissati i coefficienti del-

la (3.39), si hanno quindi n + m + 1 gradi di libertà, che possono esserefissati ad esempio minimizzando l’errore quadratico medio

E =N−1

∑n=0

[hα(n)− h(n)]2 ,


dei coefficienti della risposta impulsiva. Il problema è che in questomodo si ottiene un’equazione non lineare.

Riscriviamo la (3.39) come H(z)A(z) = B(z) ed assumiamo chehα(n) sia dato dalla risposta impulsiva di H(z−1), possiamo riscriverele equazioni nel dominio del tempo come:

hα(n) +l

∑i=1

aihα(n− i) =

{bn, per n = 0, 1, 2, . . . , m,0, per n > m,

(3.40)

ottenendo un sistema di equazioni lineari, che possono essere risolte inmodi differenti.

Approssimazioni di Padé

Nel metodo di Padé, si usano le equazioni per n > m per risolvererispetto ai coefficienti ak. Successivamente vengono risolte le prime mequazioni in cui i coefficienti ak sono adesso noti.

Approssimazioni di Prony

Nel metodo di Prony i coefficienti ak si possono trovare dalla souzionedelle equazioni normali della minimizzazione dell’errore quadraticomedio. I coefficienti bk, vengono poi determinati con lo stesso metododi Padé.

Approssimazioni di Shank

Invece di forzare la funzione a passare per i primi n + 1 punti della ri-sposta impulsiva, si considera un’interpolazione che minimizzi l’errorequadratico medio. In questo modo ricavo ancora una volta i coefficien-ti ak, dai quali si ottengono poi i coefficienti bk richiedendo che l’erroresia nullo (interpolazione passante).

4 A P P L I C A Z I O N I E T E C N O LO G I E

I sistemi frazionari possono essere realizzati con componenti frazio-nari, oppure con uno dei metodi di approssimazione descritti. Se leapprossimazioni risultano spesso il metodo più rapido per implemen-tare le nuove teorie, molteplici vantaggi si ottengono dall’implementa-zione diretta dei componenti. Il problema può sembrare complesso, sicerca di modellare con dei componenti concentrati dei sistemi che perloro natura sono diffusi, ma in realtà si è rivelato piuttosto semplice,almeno in linea di principio.

Dopo aver visto un cospicuo assortimento di basi teoriche, passiamoa qualche esempio applicativo. In questo caso, ciascun esempio spessonecessita di uno studio particolare, per comprenderne la descrizionefisica e le variabili in gioco, per questo gli esempi non saranno molti,ma saranno visti soprattutto gli esempi più noti ed immediati.

4.1 ELEMENTI FRAZIONARI CONCENTRATILa ricerca per la realizzazione di elementi elettronici a costanti concen-trate di ordine frazionario, cioè le frattanze, parte spesso dalla fisica odalla chimica. Attualmente questo tipo di componenti non è molto dif-fuso, ma potrebbero diventarlo: lo sviluppo di tali componenti infattiè di interesse economico per la registrazione di brevetti.

Le tecnologie attuali non permettono di fare un componente concen-trato che pur basandosi sullo stesso principio, sia in grado di compor-tarsi sia da induttore che da condensatore frazionario, ma, ovviamente,i due casi vanno studiati separatamente.

4.1.1 Condensatore frazionario

Pochissimi articoli trattano l’argomento, partiamo da [8] per introdurrel’argomento. Un condensatore frazionario, dovrà avere una risposta infrequenza del tipo

1sαC

, ⇔ 1(jω)αC

,

con α ∈ (0, 1).Ottenere una risposta di questo tipo è più semplice di quanto non I condensatori ideali

non esistono!possa sembrare, per assurdo, forse è proprio più complesso ottenereα = 1 che non un valore α < 1. Infatti un condensatore ideale, con undielettrico ideale non può esistere, poiché sarebbe non causale.

Un buon condensatore, uno qualunque tra i tanti disponibili in com-mercio, ha una risposta con α che si avvicina all’unità ma non è maipari ad 1. Dei valori realistici potrebbero essere da 0.999 o 0.9999, inogni caso dei condensatori reali mostrano sempre un certo grado difrattanza.

65

66 APPLICAZIONI E TECNOLOGIE

Figura 7: Parte reale della costante dielettrica in alta frequenza, in funzionedella temperatura. Si osservi l’ampiezza del range di misura, cheoscilla tra i −200 °C ed i 200 °C. Il grafico è preso da [37]

Invece di spingere gli studi, verso la ricerca dell’isolante che piùUn dielettricomodesto permette di

realizzarecondensatori

frazionari,

si avvicini alla idealità, occorre cercare materiali che si discostino diproposito dal comportamento ideale di isolante perfetto, in questomodo è possibile ottenere funzioni di trasferimento con vari gradi difrazionalità.

Una delle maggiori difficoltà residue è quella di trovare un materialecon un comportamento uniforme su un’ampia banda di frequenze etemperature.

Il Solfato Idroazotato di Litio (LiN2H5SO4) [37] è uno dei compostiun buon materiale èil solfato idroazotato

di litio!chimici che presenta le caratteristiche desiderate1, con costante dielet-trica ε(jω) = εr + jεi sia reale (εr) che immaginaria (εi) molto ampie(> 10× 106) e costanti su un ampio intervallo, come mostrato nellafigura 7, e con un valore di α ≈ −1/2.

Calcolando l’impedenza di un condensatore formato da due faccepiane parallele di superficie S e distanzate d,

Z =d

jωε0ε(jω)S,

in cui d è la distanza tra le facce, S è la superficie delle armature e lacostante dielettrica

ε(jω) = εrω−12 − jεiω

− 12 ,

per εr = εi = ε, servendosi della notazione polare dei numeri comples-si e ricordando che εj π

2 = j,

ε(jω) = εω−12 (1− j)

= εω−12√

2e−j π4 =√

2εω−12(ej π

2)− 1

2 =√

2ε(jω)−12 ,

si ottiene l’equazione

Z =d

Sjωε0ε(jω)√

2(jω)−12

=d

Sε0ε√

2(jω)+ 12

=1

CF(jω)12

(4.1)

1 Cioè presenta un comportamento pressoché costante su un’ampia banda di temperaturee frequenze. Il materiale deve essere usato come dielettrico di un condensatore, chequindi avrà un comportamento frazionario, sarà cioè una frattanza.

4.1 ELEMENTI FRAZIONARI CONCENTRATI 67

Strato dielettrico

Strato resistivo Equivalente circuitale

a b

c

a b

c

Figura 8: Schema realizzativo della capacità frazionaria a film sottile e schemaelettrico equivalente

che esprime l’impedenza (o frattanza) di un condensatore di ordinefrazionario in cui la sua capacità di ordine frazionario CF, sarà la partecostante della (4.1):

CF =√

2Sε0ε

d.

Un solo elemento così realizzato, può sostituire, e con grande accura-tezza, una delle lunghe reti a scala sintetizzabili dalle approssimazionidate al § 2.4 nella pagina 27!

Lo studio di materiali come il Solfato Idroazotato di Litio o suoiequivalenti, è ancora agli inizi, probabilmente perché un materiale delgenere in passato sarebbe stato scartato dal punto di vista elettroni-co, perché incapace di realizzare un buon isolante oppure un buonconduttore.

Condensatore a film sottile

Roy [35], propone vari metodi per la realizzazione di un impedenzafrazionaria, tra i quali una serie di reti a scala, ed un condensatore afilm sottile.

Il condensatore a film sottile realizza di fatto una rete a scala RC Reti a cascata RCpossono essererealizzate con uncondensatore a filmsottile

in cui i componenti non sono discreti ma diffusi. Lo strato resistivorealizza la serie di resistenze, lo strato di dielettrico, sotto il quale èposto uno strato conduttivo si comporta da condensatore.

4.1.2 Induttore frazionario

Sul lavoro fatto per la realizzazione di un induttore frazionario, si sasolamente che nel novembre 2006 è stato presentato un brevetto alriguardo. Il brevetto citato da [16, pagina 87] è US 20060-267597.

Si suppone, che sia possibile realizzare un induttore frazionario im-mergendo un avvolgimento elicoidale in un mezzo parzialmente iso-lante, in un modo molto simile a quanto fatto per il condensatore fra-zionario: si sostituisce il mezzo con proprietà isolanti che si avvicinanoal caso ideale, con un mezzo in grado di realizzare una funzione di tra-sferimento sα, con α ∈ (0, 1) in un ampio intervallo di frequenze etemperature.


R0

R4R6 R2

C1C3C5

Figura 9: Rete a scala RC, come modello di una linea di trasmissione o comestruttura a componenti discreti a costanti concentrate.

4.2 LINEE DI TRASMISSIONEUna sequenza di resistenze serie e condensatori (normali, di ordineLe linee di

trasmissione sono deicircuiti

“naturalmente”frazionari

intero) è naturalmente descritto da un modello frazionario [27]. La se-quenza di elementi circuitali infatti, si può scrivere facilmente con unosviluppo in frazioni continue, che, come abbiamo visto, corrispondead un sistema non intero.

È per altro noto dai corsi di campi elettromagnetici che una linea ditrasmissione può essere modellata con una rete a scala di elementi RCconcentrati, sebbene sia per sua natura un sistema distribuito. Cometutti i sistemi distribuiti trova una descrizione accurata e naturale con idifferintegrali non interi, vediamo come, seguendo il filo logico di [27].

Consideriamo il semplice circuito in figura 9, partendo dai solo dueelementi di chiusure del circuito, cioè R0 e C1. L’ammettenza di ingres-so Y1 sarà data dal parallelo, cioè dalla somma delle ammettenze:

Y1 = G0 + sC1.

Aggiungendo la resistenza R2, ho un’impedenza totale di ingressodei tre componenti

Z2 = R2 + Z1 = R2 +1

Y1= R2 +

1

sC1 + 1R0

,

che scriviamo dividendo ambo i membri per R2, e dividendo numera-tore e denominatore della frazione per C1, quindi

Z2

R2= 1 +

1R2C1

s + 1C1R0

. (4.2)

Proseguendo con la stessa semplicissima logica, il circuito formatocon quattro componenti, R0, C1, R2 e C3, avrà un’ammettenza d’ingres-so:

Y3 = sC3 + Y2,

ed aggiungendo R4 si ha poi

Z4 = R4 +1

Y3= R4 +

1

sC3 + 1Z2

.

Come per l’impedenza Z2 che è stata riscritta nella forma (4.2), riscri-viamo Z4 nel modo seguente:

Z4

R4= 1 +

1R4C3

s + 1C3Z2

, (4.3)

4.2 LINEE DI TRASMISSIONE 69

poi dividiamo per R2 l’ultimo quoziente parziale

1C3Z2

=1

R2

C3Z2R2

,

in modo da poter sostituire z2/R2 con l’espressione esplicita (4.2).La (4.3) in definitiva si scrive

Z4

R4= 1 +

1R4C3

s +1

C3R2

1 +1

R2C1

s + 1C1R0

, (4.4)

espressione perfettamente analoga alla (4.2).Definendo le frequenze

ω0 =1

C1R0ω1 =

1R2C1

ω2 =1

C3R2ω3 =

1R4C3

,

possiamo riscrivere l’impedenza di ingresso (4.4) per la rete a scala dicinque componenti nella forma più compatta

Z4

R4= 1 +

ω3

s+ω2

1+ω1

s+ω0,

grazie anche alla notazione sintetica di Rogers delle frazioni continue2.La serie ottenuta, che può essere facilmente estesa allo stesso modo

per un numero maggiore di elementi, inoltre si dimostra essere con-vergente [44]. Per ottenere una forma standard di sviluppo in frazio-ni continue facciamo un’altro passo, moltiplichiamo ciascuna frazioneparziale con denominatore s+, per il termine 1/s, nella forma estesa siottiene

Z4

R4= 1 +

ω3

s+ω2

1+ω1

s+ω0 = 1 +

ω3s

1 +ω2s

1 + ω11s

(s+ω0) 1s

,

poi, con la sostituzione

ωis

= ai,

otteniamo

Z4

R4= 1 +

a3

1 +a2

1 + a11+a0

,

cioè una frazione continua in forma canonica.Per completezza scriviamo la generalizzazione al caso di N compo-

nenti, con N pari, si avrà

ZNRN

= 1 +aN−1

1+aN−2

1+· · · a2

1+a1

1+a0 (4.5)

2 Vedi § 2.4.1 nella pagina 28 per maggiori informazioni.


R0

RN Ri-1

C1CiZN

i=3,5,7...N-1

Figura 10: Schema generico di una rete a scala, con N componenti. Il bloccoevidenziato in azzurro è ripetuto più volte fino ad avere in totale ilnumero di componenti desiderati

in cui si è fatto riferimento alla figura 10.La somma di tutti i termini CFE, per

R0 = R2 = R4 = · · · = RN−2 = R, RN =12

R,

C1 = C3 = C5 = · · · = CN−1 = C,

possiede uno sviluppo convergente, come mostrato da [44], pari a

a1+

a1+· · · a

1+a

1+a =

√4a + 1

1 +(√

4a+1−1√4a+1+1

)N + 1

−√

4a + 12

− 12

.

Combinando lo sviluppo dell’espansione CFE con la (4.5), e tenendoconto dei valori (4.2), si ottiene infine

Z

√Cs

R=

√4a + 1

4a

[(√4a + 1 + 1

)N+1(√4a + 1− 1)N+1(√

4a + 1 + 1)N+1(√4a + 1− 1

)N+1

],

cioè un’espressione dell’impedenza di ingresso della linea (o della reteL’impedenza èproporzionale alla

radice della frequenzaa scala RC), proporzionale alla radice del termine a, che rappresental’inverso del prodotto tra resistenza e capacità. La natura fraziona-ria del sistema discende proprio dalla dipendenza dalla radice dellafrequenza.

4.3 PRESTAZIONI ALGORITMI NUMERICIProviamo a vedere il comportamento delle approssimazioni mostrateStudio di un sistema

test al capitolo 3 nella pagina 39, su un sistema test, ad esempio possiamoconsiderare

y(t) = 0.7 aD0.3t x(t) + 0.4 aD

0.7t x(t), (4.6)

ottenuto dall’equazione (2.1) ponendo:

α = {α0, α1} = {0.3, 0.7},N = 0, M = 1.

Scegliamo un segnale di forma triangolare, cioè x(t) = tri(t), con

−1 0 10

1tri(t) ampiezza unitaria, origine in a = −1 e fine in b = 1, che potremmo

campionare con periodo

T =1

100


L’analisi dei differenti metodi è stata eseguita in Matlab, per cia-scuno dei metodi di approssimazione proposto è stato realizzato unoscript ad hoc per il tracciamento delle differenti curve di approssima-zione. Per compattezza negli script si è studiata solo la prima rampacrescente del segnale per la derivata di ordine 0.3, le considerazioni egli script proposti permettono di estendere facilmente l’applicazione asistemi più complessi. L’inizializzazione è mostrata nel listato 4.1.

1 % Esempio di simulazioni per il confronto di metodi numerici.

2

3 % Studio il sistema y(t) = D0.3 {x(t)}, con ingresso rampa tra

[-1,0]

4 a=0.3; % ordine di derivazione \alpha

5 T=1/100; % periodo

6 t=-1:T:0; % dominio del tempo

7 N=50; % numero di punti usati dalle approssimazioni

8 x=t+1;Listato 4.1: Script per il confronto dei differenti algoritmi di elaborazione diordine non intero. Inizializzazione dei parametri.

4.3.1 Soluzione analitica

La soluzione analitica si può ottenere dalla definizione di Riemann- La soluzioneanalitica si puòottenere risolvendol’integrale di R-L

Liouville (1.14) nella pagina 7, che quindi si scriverà

−1D0.30 x(t) =

1Γ(−0.3)

∫ t

−1(t− τ)−0.3−1(τ + 1)dτ

nel dominio t ∈ [−1, 0]. Spezziamo poi l’integrale nella somma di dueintegrali∫ t

τ=−1(t− τ)−1.3τdτ +

∫ t

τ=−1(t− τ)−1.3dτ per t ∈ [−1, 0]

che risolveremo separatamente, ricordando poi di reintrodurre allafine il termine 1/Γ(−0.3) tralasciato.

Primo integrale

Con le sostituzioni

t− τ → u ⇒ τ = t− u,

dτ = −du,

per τ = −1, u = t + 1,

per τ = t, u = 0,

possiamo risolvere il primo integrale∫ t

τ=−1(t− τ)−1.3τdτ = −

∫ 0

u=t+1u−1.3(t− u)du

= −t∫ 0

u=t+1u−1.3du +

∫ 0

u=t+1u−1.3u du

= t[

u−0.3

0.3

]0

u=t+1+[

u0.7

0.7

]0

u=t+1

= −t(t + 1)−0.3

0.3− (t + 1)0.7

0.7(4.7)


con i semplici passaggi matematici illustrati.

Secondo integrale

Analogamente per il secondo integrale, con le stesse sostituzioni, pos-siamo scrivere∫ t

τ=−1(t− τ)−1.3dτ = −

∫ 0

u=t+1u−1.3du = − (t + 1)−0.3

0.3, (4.8)

con passaggi identici al caso precedente.

Ricomposizione della soluzione

Rimettiamo quindi insieme tutti i pezzi, cioè la somma delle soluzionidei due integrali (4.7), (4.8) moltiplicati per la funzione gamma, e ricor-dando inoltre che il dominio della variabile t è compreso nell’intervallo[−1, 0], per cui

−1D0.30 x(t) = − 1

Γ(−0.3)

[(t + 1)0.7

0.7+ t

(t + 1)−0.3

0.3+

(t + 1)−0.3

0.3

]= − 1

Γ(−0.3)

[(t + 1)0.7

0.7+

(t + 1)−0.3

0.3(t + 1)

]otteniamo infine l’equazione analitica del sistema

−1D0.30 x(t) = − (t + 1)0.7

Γ(−0.3)

(1

0.3+

10.7

)(4.9)

Il calcolo della (4.9) si ottiene banalmente dalle seguenti istruzioni

12 % SOLUZIONE ESATTA, calcolata dalla definizione di Riemann

Liouville

13 % L’uscita è data dalla seguente formula

14 y1 = -(t+1).^0.7/gamma(-a) * (1/0.3 + 1/0.7);

15 plot(t,y1,’LineWidth’,2);

16 hold all

17 gridListato 4.2: Script per il confronto dei differenti algoritmi di elaborazione diordine non intero. Soluzione esatta del sistema.

l’altra parte della soluzione si può ottenere semplicemente sostituen-do gli esponenti e mettendo in evidenza un segno negativo, per questonon sarà qui riportata.

4.3.2 Metodi alle differenze

Differenze di Grünwald-Letnikov

Vediamo ora come affrontare praticamente lo stesso sistema (4.6) conil metodo proposto nel § 3.4 nella pagina 49, basato sulla definizionedi differintegrale di Grünwald-Letnikov (1.13).

Scriviamo l’equazione del sistema (4.6), in base all’espressione (3.25),ottenendo

y[n] = 0.7 ΛT0.3 · x + 0.4 ΛT

0.7 · x


in forma vettoriale, oppure

y[n] = 0.7L

∑k=0

Λ0.3[k] x[n− k] + 0.4L

∑k=0

Λ0.7[k] x[n− k]

in forma estesa, con gli elementi esplicitati rispetto all’indice k, in cuiinoltre

Λ0.3[k] =Γ(k− 0.3)

Γ(−0.3)Γ(k + 1)Λ0.7[k] =

Γ(k− 0.7)Γ(−0.7)Γ(k + 1)

,

come detto nel §3.4.In questo caso sono state realizzate due funzioni per eseguire l’algo-

ritmo: la prima, vg, calcola i coefficienti del vettore gamma Λ, è stataintrodotta nel listato 3.1 nella pagina 52; la seconda, df, esegue la con-voluzione su una finestra del segnale, con i coefficienti della funzionegamma. Vediamo quindi per prima cosa la funzione df per le derivatefrazionarie:

1 function y=df(x,q,N)

2 % Controllo ingressi.

3 if (nargin ~=2) & (nargin ~=3)

4 error(’Chiamata con numero di argomenti non valido’)

5 end

6

7 if length(x) < N

8 error(’Il vettore x deve contenere più punti degli N usati’)

9 end

10

11 if nargin==2

12 N=100;

13 end

14

15 % Inizializzazione

16 G=vg(q,N);

17 X=zeros(1,N);

18 % X(1)=x(1);

19 % y=x(1);

20

21 for i=1:length(x)

22 X(2:N) = X(1:N-1);

23 X(1) = x(i);

24 y(i) = G * X’;

25 endListato 4.3: Script per il calcolo delle derivate di ordine non intero a partiredal vettore Λ dei coefficienti binomiali ottenuti con la funzionegamma.

Lo script per calcolare l’uscita della funzione considerata con l’ap-prossimazione delle derivate con il metodo illustrato, si riduce quindialle istruzioni:

20 % DERIVATE CON LA DEFINIZIONE DI GRUNWALD, approssimata per N punti

21 % Calcolo l’uscita con la convoluzione del vettore gamma

22 y2=df(x,0.3,N) * T^(-a);

23 plot(t,y2);Listato 4.4: Script per il confronto dei differenti algoritmi di elaborazione diordine non intero. Metodo dell’approssimazione delle derivate conla definizione di G-L.


26 % METODO DI EULERO, l’equazione alle differenze in z^-1 è

27 % approssimata; si ottengono i coefficienti della definizione di G-

L

28 hE=vg(a,N);

29 y3=conv(hE,x);

30 y3=y3(1:length(t))* T^(-a); %tronco y ai soli valori della funzione

31 plot(t,y3,’LineStyle’,’--’,’color’,color(5))Listato 4.5: Script per il confronto dei differenti algoritmi di elaborazione diordine non intero. Metodo di Eulero.

Metodo di Eulero

Il metodo di Eulero parte da presupposti differenti rispetto al metodoIl metodo di Eulero èdiverso nelmetodo. . .

visto, ma poi giunge agli stessi risultati. Le basi del metodo, ed ilconfronto con l’approssimazione delle derivate con la definizione diG-L è mostrata nel § 3.5.3 nella pagina 57.

Procedendo per gradi per prima cosa si trasforma il sistema scriven-do

Y(s) = 0.7s0.3X(s) + 0.4s0.7X(s),

poi si discretizza la variabile s di Laplace con la (3.30), e si sviluppa ilbinomio risultante(

1− z−1)α

caratteristico del metodo, ottenendo

Y(z) = 0.7H0.3(z−1)X(z) + 0.4H0.7(z−1),

la quale, infine, può essere invertita calcolandone la Z −1 per tornareall’equazione in n

y[n] = 0.7h0.3[n]⊗ x[n] + 0.4h0.7[n]⊗ x[n]

che riporta le derivate ad una convoluzione tra il segnale n[n] eduna risposta impulsiva hα[n] che dipende dall’ordine di derivazioneα considerato.

La risposta impulsiva hα[n] coincide con il vettore gamma Λα intro-. . . ma giunge aglistessi risultati! dotto al § 3.4 nella pagina 49, per questo può essere calcolata nello

stesso modo, facendo uso della funzione lsvg.In questo caso la convoluzione è stata lasciata alla funzione conv

presente in Matlab, il codice necessario è mostrato nel listato 4.5.

Metodo di Tustin

Consideriamo l’approssimazione trapezoidale studiata nel § 3.5.4 nel-la pagina 59, per ottenere i coefficienti della risposta impulsiva (3.35)del sistema. Il procedimento formale è lo stesso del metodo di Euleroappena visto, l’unica differenza tangibile risiede nei coefficienti del-la risposta impulsiva (3.35), che sono ottenuti a partire da un’appros-simazione trapezoidale (3.34), invece che dalle differenze all’indietrocaratteristiche del metodo di Eulero.

La risposta impulsiva (3.35) qui riportata

hα[n] =(

2T

)α n

∑k=0

(−1)k(

α

n

)(−α

n− k

), (3.35)


34 % METODO DI TUSTIN

35 h1=vg(a,N); % La prima parte si calcola allo stesso modo

36 hT=zeros(1,N);

37 for n=1:N,

38 h2=ones(1,n); % secondo coefficiente binomiale

39 for i=1:n-1,

40 h2(i+1)=h2(i) * (i-1-a)/i * (-1)^i;

41 end

42 for i=1:n, % sommatoria delle parti

43 hT(n)=hT(n) + h1(i) * h2(end-i+1);

44 end

45 end

46

47 y4=conv(hT,x); % convoluzione troncamento e scaling

48 y4=y4(1:length(t))* (1/((2+4*a)*T))^(a);

49

50 plot(t,y4,’color’,color(3))Listato 4.6: Script per il confronto dei differenti algoritmi di elaborazione diordine non intero. Metodo di Tustin.

è stata calcolata in due fasi. Come risolvere ledifficoltà del calcolodei coefficientibinomiali

1. Sono calcolati tutti i termini

(−1)k(

α

n

)della sommatoria, per k = 0, 1, . . . , N, in cui N è il numero di pun-ti considerati nell’approssimazione. Ciò si ottiene semplicementecon h1=vg(a,N).

2. È stata calcolata la seconda parte(−α

n− k

)della sommatoria, per k = 0, 1, . . . , n per ciascuno degli indicin considerati in uscita. Si noti come in questo caso l’estremosuperiore n vari, così come il numero di elementi per ciascuncampione n della risposta impulsiva. I termini sono stati calco-lati con una procedura iterativa, analoga a quella utilizzata nellafunzione vg, e mostrata nel § 3.4 nella pagina 49, e memorizzatinel vettore h2.

3. Dagli elementi parziali della sommatoria, memorizzati nelle va-riabili h1, h2, si è composta la funzione di trasferimento hT carat-teristica del metodo, con le istruzioni del 3° ciclo for tra le righe42–44.

Metodo di Al-Alaoui

Il metodo di Al-Alaoui è stato presentato nel § 3.5.5 nella pagina 62,vediamo adesso come applicarlo.

In modo analogo rispetto al procedimento applicato per il meto-do di Tustin, la sommatoria che esplicita i coefficienti della rispostaimpulsiva, è spezzata in due parti, poi sommate separatamente.


0 5 10 15 20 25 30 35−0.06

−0.05

−0.04

−0.03

−0.02

−0.01

0

Confronto h[n]

h[n]

EuleroTustinAl-Alaoui

Figura 11: Confronto tra i coefficienti della risposta impulsiva di ciascuno deimetodi studiati.

In generale, i coefficienti binomiali si calcolano tutti nel modo in-dicato al § 3.4 nella pagina 49; i tentativi di calcolare i coefficientibinomiali dalla definizione diretta (1.5), per semplice sostituzione delfattoriale con la funzione gamma, non permettono di ottenere risulta-ti accettabili, sebbene per approssimazioni con meno di 50 elementisi stia comodamente all’interno dell’insieme dei valori rappresentabilidai dati di tipo double.

Ciascun coefficienti binomiali quindi, è stato calcolato con un al-goritmo iterativo, così come fatto per gli elementi del vettore Λ. Laporzione di script che restituisce la risposta impulsiva del metodo èmostrata nel listato 4.7 risultati sono riepilogati nel codice:

L’andamento generale delle risposte impulsive dei metodi appareuniforme, in tutti i casi si ha un primo coefficiente unitario ed unarapida decrescita del modulo dei coefficienti, che diventano via viameno determinanti per l’uscita, man mano che l’indice v avanza. Unconfronto tra i diversi metodi è mostrato in figura 11.

4.3.3 Riepilogo risultati

Le curve ottenute con i vari metodi sono mostrati nella figura 12 afronte, in relazione al segnale x(t) considerato, ed in relazione allasoluzione esatta, indicata con un tratto più spesso.

Le risposte nel dominio del tempo non mostrano alcuna caratteristi-ca di rilievo, gli algoritmi si comportano tutti più o meno allo stessomodo, per apprezzare le differenze è necessario infatti passare al do-minio della frequenza, in cui, dal mappaggio tra il piano s ed il pianoz possiamo ricavare le prestazioni dei diversi algoritmi.


−1 −0.8 −0.6 −0.4 −0.2 00

0.2

0.4

0.6

0.8

1

1.2

1.4Confronto algoritmi

t [sec]

y

EsattaGrunwaldEuleroTustinAl-Alaoui

Figura 12: Confronto tra le differenti approssimazioni citate. In blu il segna-le in ingresso al sistema, con tratto spesso la soluzione esatta diriferimento per tutte le approssimazioni.

−1 −0.8 −0.6 −0.4 −0.2 0−0.08

−0.06

−0.04

−0.02

0

0.02

0.04

Errore in uscita

t [sec]

GrunwaldEuleroTustinAl-Alaoui

Figura 13: Ampiezza dell’errore dei vari metodi di approssimazione utilizzati.


54 hA1=vg(a,N); % Primo coefficiente binomiale

55 hA=zeros(1,N);

56 for n=1:N,

57 %hn=0;

58 hA2=ones(1,n); % secondo coefficiente binomiale

59 for i=1:n-1,

60 hA2(i+1)=hA2(i) * (n-i+a)/(n-i+1) * (-1)^i;

61 end

62 for i=1:n, % sommatoria delle parti

63 hA(n)=hA(n) + hA1(i) * hA2(end-i+1) * (1/7)^(n-i+1);

64 end

65 end

66 hA=hA/hA(1);

67 y5=conv(hA,x); % convoluzione troncamento e scaling

68 y5=y5(1:length(t))* (8/(7*T))^(a);

69

70 plot(t,y5,’color’,color(4))Listato 4.7: Script per il confronto dei differenti algoritmi di elaborazione diordine non intero. Metodo di Al-Alaoui.

4.4 ESEMPI APPLICATIVIConcludiamo il lavoro con degli esempi più concreti, cercando di pas-sare dall’astratto a semplici casi applicativi. Tuttavia, ciascun caso rea-le richiede uno studio mirato al problema, per cui qui si andrannoad analizzare solamente due esempi pratici, un oscillatore smorzato,ed un sistema di sintesi del suono basato sul ritardo frazionario. Sirimanda a lavori futuri il compito di elaborare e sviluppare approfon-ditamente gli argomenti proposti, o tutti gli altri solo accennati o noncitati affatto, per i quali il calcolo frazionario risulta essere un potentestrumento di sviluppo.

4.4.1 Oscillatore smorzato

Un sistema dinamico che si presta facilmente ad uno studio sistema-tico è un sistema risonante attenuato, meccanicamente costituito daUn tipico esempio di

studio: risonatoresmorzato

una massa collegata ad una molla, immersa in un fluido in grado dismorzare le vibrazioni. Le possibili varianti fisiche del sistema descrit-to sono molteplici, le equazioni possono essere adottate con facilità,a qualsiasi sistema elettronico o meccanico, che vibri o oscilli: al ca-so di corde vibranti, pendoli, sistemi in rotazione fuori asse, sistemielettronici risonanti elementari, e molti altri.

È possibile considerare un sistema non intero di questo tipo, se-guendo l’approccio di [22], in cui lo smorzamento è di ordine α nonintero:

mD2x(t) + b Dαx(t) + kx(t) = f (t).

Introducendo i parametri di frequenza naturale (frequenza di oscil-lazione del sistema non smorzato)

ω0 =

√km


Tabella 3: Radici dell’oscillatore smorzato frazionale, ricavate con il comandoroots di Matlab, per ω0 = 1 rad s−1, e ζ = 0.05.

Radice Valore

ζ1 0.7073 0.7319 iζ2 0.7073 −0.7319 iζ3 −0.7073 0.6819 iζ4 −0.7073 −0.6819 i

e smorzamento

ζ =b

2mω2−α0

,

è possibile riscrivere l’equazione del sistema nella forma

D2x(t) + 2ω2−α0 ζ Dαx(t) + ω0x(t) = f (t).

Applicando poi l’usatissima (1.26), si può ottenere la funzione di tra-sferimento

H(s) =F(s)X(s)

=1

s2 + 2ω2−α0 ζ + ω2

0,

il cui polinomio indiciale, fissato

α =12

coincidente con il massimo comun divisore degli esponenti, e riscrittoin funzione di multipli di α indicati con ς = sα

ς4 + 2ω2−α0 ζς + ω2

0 = 0,

deve essere risolto per trovare i poli, ad esempio tramite una procedurastandard.

In Matlab è possibile trovare le radici di un polinomio dalla matricedei coefficienti c con il comando

1 roots(c)

Se fissiamo i valori

ω0 = 1 rad s−1, ζ = 0.05,

la matrice dei coefficienti c sarà [1, 0, 0, 0.1, 1], il comando roots

restituisce le radici mostrate nella tabella 3.Per tornare in s, la relazione s = ς1/α inversa della trasformazione

ς = sα, raddoppia la fase delle radici calcolate. Le radici ς1 e ς2 so-no caratterizzate da una fase il cui modulo non supera i π

2 rad, percui i rispettivi poli in s con fase doppia (1/α) rispetto a ς sarannocaratterizzati da arg(s1,2) < |π|, e quindi saranno localizzati sulla su-perficie principale di Riemann. Le altre due soluzioni non rispettanoquesta condizione, non hanno senso in un sistema fisico reale e sonoda scartare.

Invertendo la funzione trovata si ottiene la risposta impulsiva delsistema che presenta sia parte reale che immaginaria. Estraendo lasola parte reale, la risposta impulsiva

h(t) = <{

r1ς1etς21 u(t) + 0Dα

t r1etς21

},

ci fornisce l’andamento temporale delle oscillazioni del sistema.


4.4.2 Codifica del parlato

Le peculiarità sonore del parlato, permettono di suddividere tutti isuoni emessi, principalmente in due categorie, i suoni con intonazione,come le vocali, ed i suoni rumorosi, come molte consonanti, quali adesempio s, t, c. . . .

Le vocali sono formate da suoni intonati, si ha quindi una regolari-Le vocali sonocostituite da onde

sonore periodiche. . .tà periodica del segnale sonoro, che in frequenza sarà caratterizzatopicchi di ampiezza in corrispondenza delle frequenze fondamentali oarmoniche3 della voce.

I sistemi impiegati per il parlato adottano frequenze di campiona-mento tipicamente di 8 KHz [23]. Se si considera una voce femminilecon frequenza fondamentale di 250 Hz quindi con periodo 1/250 =4 ms, il numero di campioni per periodo sarà pari a 4 · 8 = 32 samples.Se il sistema di formazione e predizione dovesse introdurre un erroredi un campione, si avrebbe un’alterazione di tono circa del 3 % (1/32),ben percepibile ad orecchio, con conseguente perdità di qualità ed ac-curatezza del sistema. Tra l’altro, la frequenza fondamentale indicatadi 250 Hz non è affatto acuta, per frequenze maggiori l’errore cresce!

I sistemi per la formazione dei suoni sono illustrati ad esempio negliarticoli [19, 18], una sintesi è disponibile su [23], in cui è invertitol’ordine dei sottosistemi di predizione del tono e generazione delle formanti.Questa sezione non è intesa come parte di approfondimento alla sintesidel parlato, l’argomento è sufficientemente vasto da costituire un temaper altre tesi e tesine, si vuole solamente illustrare in linea generaleun modo possibile per formulare un sistema alle differenze non intere.Un’esempio applicativo in Matlab si può trovare in [38].

Sempre concentrandoci sulle vocali (sonore), cerchiamo un modelloche simuli in qualche modo i risonatori del tratto vocale, e la sorgenteacustica responsabile dell’eccitazione impressa. Costruiamo il sistemain due parti separate, una sorgente opportuna ed un modello del trattovocale risonante.

Il modello deve essere a coefficienti dinamici, poiché durante la fo-. . . l’andamento deltono rende il senso

del discorso connaturalezza. . .

nazione i diversi suoni con i quali sono costruite le parole ovviamentecambiano, per via dell’azione muscolare sui diversi apparati coinvoltinella costruzione delle parole. I cambiamenti sono lenti rispetto adun sistema elettronico, questo permette di aggiornare i coefficienti deifiltri (filtri adattativi) ad intervalli di circa 20 ms.

Basandosi sui modelli di [19, 18] possiamo porre il pitch predictor perprimo, come mostrato in figura 14 nella pagina successiva. Lo stimato-re del tono prende in ingresso un segnale simile ad un ronzio, che nelcaso di una vocale può essere composto da una sequenza di impulsi,generati da un’altra parte del circuito che occupa di effettuare l’anali-si del parlato sintetizzato, in un’anello di controreazione. Il filtro disintesi a sua volta è costituito da un pitch predictor, che si occupa dellastima di fino del tono, consentendo la modulazione del discorso, e daun formant filter che modella i suoni (shaping in frequenza) conferendole caratteristiche timbriche e fonetiche desiderate al discorso.

Tipicamente il formant filter è rappresentato in frequenza da una ri-sposta con 3 o più picchi, ciascuno dei quali corrisponde ad una dellearmoniche principali, la distribuzione e l’ampiezza delle armonichedetermina il timbro esatto, quindi determina la vocale sintetizzata.

3 Chiaramente, le armoniche caratterizzano il timbro della voce, e saranno situate amultipli interi della frequenza fondamentale che determina l’intonazione.


Pitchpredictor

Formantfilter

S Filtro di sintesi

Parlato di qualità

Sorgente non modellata

Figura 14: Schema di un sistema di sintesi del parlato. Il segnale di ingressodi natura rumorosa e non modellato, simile ad un ronzio, passa nelfiltro di predizione di pitch che stabilisce il tono con precisione, poiviene modellato dal filtro formante.

Il pitch predictor deve modulare il tono secondo l’andamento della . . . un filtro nonintero permette undiscorso più fluido!

parola o del contesto della frase in cui si trova, per rendere fluido escorrevole il discorso. Come accennato all’inizio di questa sezione, se iltono fondamentale è già acuto, è impossibile ottenere transizioni fluidee non evidenti tra i diversi toni che si susseguono nella parola, ma siotterrebbe invece un’andamento a scala. L’uscita del pitch predictorè una sequenza di impulsi, per ottenere transizioni tonali più fluide ènecessario variare con continuità la spaziatura (periodo) tra gli impulsi,per questo è possibile usare un sistema alle differenze non intere.

Lo studio di sistemi di sintesi del parlato di ordine frazionario sista diffondendo, alcuni studi effettuati dimostrano come sia possibi-le ottenere un miglioramento del Rapporto Segnale Rumore (SNR) eduna riduzione del numero di bit necessari. Per una bibliografia diriferimento sul tema si rimanda a [23].

4.4.3 Conclusioni

La difficoltà primaria di questo lavoro è stata la ricerca e lo studiosistematico di tutte le teorie alla base del calcolo frazionario.

Il duro impatto iniziale per apprendere l’argomento, servendosi de-gli articoli specializzati disponibili, è stato lentamente superato, soltan-to con il tempo e la perseveranza: ciascun articolo infatti si presentavaparzialmente incompleto, dando per note le nozioni di base necessa-rie alla comprensione del testo, e appoggiandosi pesantemente ad unvasto materiale bibliografico dello stesso genere; non solo: le infor-mazioni fornite, sempre estremamente sintetiche, non sono ovvie, marichiedono una successiva espansione per poter essere comprese inmodo chiaro.

Quello che si è cercato di fare, è stato proprio un lavoro di compren-sione, di ricerca e di studio, raccogliendo ed organizzando nel modopiù coerente possibile la moltitudine di informazioni preziose, ma didifficile accesso, sparse tra i vari articoli.

Per ciascun argomento trattato, sono riportate le indicazioni biblio-


grafiche della fonte, le informazioni reperite sono sempre state integra-te con considerazioni personali, i procedimenti matematici sono statiampliati per facilitare la comprensione di tutti i passaggi, sono stateproposte porzioni di codice Matlab per passare rapidamente ad unafase applicativa e vari esempi chiarificanti ideati dall’autore. Oltre agliarticoli riportati in bibliografia ne sono stati consultati numerosi altri,che si sono rivelati però infruttuosi. In bibliografia infatti sono presen-ti soltanto le fonti dalle quali sono state attinte informazioni riportatein questo testo, ma per chiunque voglia approfondire, una bibliografiamolto numerosa è disponibile a chiusura degli articoli [30, 23], o dellibro [16].

Non è stato possibile entrare nel merito di uno studio approfonditodi qualche problema specifico, si è accennato ad esempio al caso dellasintesi del parlato, o al problema degli oscillatori, ma una realizzazionefunzionante di un sintetizzatore vocale ad esempio, avrebbe richiestoun approfondimento abbastanza vasto da risultare sufficiente per unaltro lavoro di tesi.

Senza prima una buona comprensione di base dell’argomento, nonsarebbe stato possibile affrontare applicazioni specifiche in dettaglio,per cui si spera con questo testo di fornire una rapida guida introdutti-va allo studio dei sistemi frazionari, per permettere a chiunque vogliacimentarsi in una tesi sull’argomento di apprendere rapidamente lebasi e la visione d’insieme della materia.

A TO O L B O X N I N T E G E R

In Matlab sono disponibili potenti strumenti per lo studio di siste-mi elettronici di vario genere. Questi strumenti, raccolti in apposi-ti toolbox lavorano tipicamente su oggetti1 di tipo filtri o funzioni ditrasferimento, purché siano di tipo intero.

Con la venuta dei sistemi di ordine non intero, stanno comparen-do pian piano vari strumenti per poter espandere le funzionalità diMatlab alla nuova classe di sistemi, ed in particolare sono noti duetoolbox creati per lo scopo:

NINTEGER. Il toolbox ninteger è un pacchetto liberamente accessibile cheimplementa per Matlab varie funzionalità relative ai sistemi diordine non intero.

CRONE Il toolbox crone è stato sviluppato dal team di Oustaloup et al.,ma non è disponibile gratuitamente, e si basa sulle metodologiedi controllo Contrôle Robuste d’Ordre Non Entier (CRONE).

Il pacchetto crone non è disponibile gratuitamente, sebbene appa-rentemente sia più avanzato, con un’implementazione che segue ilmetodo ad oggetti adottano nel toolbox esistenti di Matlab.

Per questo motivo, daremo una rapida occhiata a ninteger, menoavanzato e attualmente basato su funzioni specifiche piuttosto che suoggetti, ma liberamente disponibile sul sito http://www.mathworks.com/matlabcentral/fileexchange/8312. Ulteriori informazioni sonodisponibili su un articolo descrittivo [42], o sul manuale d’uso delpacchetto [41].

Un ulteriore set di script è disponibile sul sito http://www.acoustics.hut.fi/software/fdtools/, si tratta di una raccolta di funzioni persistemi frazionali tempo-discreti, con il nome di fdtools.

A.1 PROGRAMMAZIONE AD OGGETTICome già accennato nell’introduzione, in Matlab sono presenti varitoolbox che si propongono come un valido sussidio per l’ingegnereelettronico, in particolare si ha:

• filterdesign, per la progettazione e l’implementazione di filtrianalogici, o digitali deterministici ed adattativi;

• controlsystem, per la progettazione e l’implementazione di siste-mi di controllo, sia analogici che digitali.

I toolbox esistenti sono basati sul metodo di programmazione ad og-getti, che permette di costruire filtri o funzioni di trasferimento, come

1 Si intendono gli oggetti della programmazione, dati strutturati di vario tipo, ciascuno deiquali possiede metodi specifici per l’implementazione delle operazioni, come vedremopiù avanti.

83



http://www.acoustics.hut.fi/software/fdtools/

http://www.acoustics.hut.fi/software/fdtools/

84 TOOLBOX NINTEGER

oggetti, ossia come classe di dati strutturati che contengono tutte le in-formazioni necessarie. Per ciascuna classe esiste uno specifico metodoche implementa in modo opportuno, dipendente dalla classe, alcuneelaborazioni per l’oggetto considerato. Nel gergo della programma-zione ad oggetti si dice che si ha un overloading dell’operatore, cioè sihanno diversi metodi selezionati dal tipo di classe con la quale si staoperando.

Chiariamo con un esempio. Supponiamo di voler conoscere le di-mensioni di un sistema. Per fare questa operazione esiste il comandosize, che può operare sia su una normale matrice, restituendo il nume-ro di elementi lungo ciascuna dimensione (righe, colonne. . . ), sia suun sistema intero, restituendo l’ordine del sistema stesso, cioè il gradomassimo.

Un sistema tra l’altro, può essere specificato in vari modi, cioè co-me funzione di trasferimento in forma di rapporto di polinomi, comeproduttoria di zeri e poli e guadagno, a ciascuna specifica i dati sonomemorizzati in modo diverso, in una struttura di dati oggetto di unparticolare tipo.

Se il sistema è creato nella forma di funzione di trasferimento conil comando tf, l’esecuzione del comando size sull’oggetto tf vieneeseguito con il metodo specifico per la classe tf. Un’altro script, conun diverso algoritmo, provvederà a restituire il risultato dell’istruzionesize su una funzione di trasferimento di tipo diverso, come ad esempioun oggetto poli zeri e guadagno ottenuto con il comando zpk.

In Matlab, il meccanismo di selezione del metodo appropriato, inbase alla classe dell’argomento della funzione, avviene in modo au-tomatico, rispettando alcune semplici regole di strutturazione delle di-rectory contenenti le funzioni: gli operatori per una data classe devonotrovarsi all’interno di una cartella con il nome formato dal carattere @seguito dal nome della classe (@nomeclasse).

Esiste quindi una funzione size, sia nella cartella control/@zpk, chenella cartella control/@tf. Per l’utente la distinzione non è palese,l’esecuzione del comando size infatti è identica in entrambi i casi!

Le informazioni fornite in questo testo sono utili per capire comeestrapolare i dati necessari dai sistemi comunemente costruiti in Ma-tlab. Gli script forniti in questo testo non vanno ad espandere colmetodo della programmazione ad oggetti i pacchetti forniti con Ma-tlab, ad ogni modo ci è sembrato opportuno dare dei cenni di baseper capire da subito come operare con gli oggetti esistenti e come poterespandere i toolbox esistenti in maniera organica.

A.2 NINTEGER: COMANDI COMUNIVediamo una breve descrizione dei comandi messi a disposizione dalpacchetto ninteger, per capire da subito come operare in pratica con isistemi frazionari. Maggiori informazioni ovviamente sono disponibilisulla guida del toolbox [41, 42]

bodeFr

Permette di generare i diagrammi di Bode per i circuiti di ordine nonintero. La sintassi richiede di cambiare il nome del comando, proprio


per l’assenza di implementazione di un oggetto di sistema di ordinenon intero. La sintassi è:

1 [gain, phase, w] = bodeFr(F, Q, w, delay)

in cui i parametri sono:

F sistema in esame. Può essere nella forma tf oppure può contenerei parametri di un sistema PID in una matrice della forma [P I

lambda D mu];

Q è l’esponente globale frazionario, cioè il massimo comun divisoretra gli esponenti, o il minimo comune multiplo;

W definisce in modo esplicito il vettore opzionale delle frequenze diinteresse;

DELAY specifica un ritardo opzionale da aggiungere al filtro.

Vedremo tra breve, come specificare un filtro di ordine frazionario

Modelli matematici dei sistemi

Se costruiamo un semplice filtro, ad esempio con le istruzioni

1 F = tf(1,[1 1])

otteniamo un oggetto

F =1

s + 1

di tipo transfer function, cioè un oggetto funzione di trasferimento. Perottenere un elenco delle proprietà specifiche di un oggetto, in questocaso della funzione di trasferimento F, basta scrivere

1 get(F)

per ottenere il seguente output

1 num: {[0 1]}

2 den: {[1 1]}

3 ioDelay: 0

4 Variable: ’s’

5 Ts: 0

6 InputDelay: 0

7 OutputDelay: 0

8 InputName: {’’}

9 OutputName: {’’}

10 InputGroup: [1x1 struct]

11 OutputGroup: [1x1 struct]

12 Name: ’’

13 Notes: {}

14 UserData: []

che mostra tutte le proprietà relative.Si osservi che in Matlab, un singolo oggetto è in grado di specifi-

care un sistema Multiple Input Multiple Output (MIMO), per cui perimmagazzinare dati irregolari in oggetti di tipo tf, sono utilizzati degliarray di tipo cell, che possono essere immaginati come particolari ma-trici in cui gli elementi possono essere dati difformi di qualsiasi tipo.

86 TOOLBOX NINTEGER

Per accedere agli elementi di un dato di tipo cell, si usano le parentesigraffe, mentre per accedere alle singole proprietà di un dato struttu-rato, si usa il punto seguito dal nome della proprietà. In definitivaquindi, i coefficienti del denominatore del sistema monodimensionaleF precedentemente indicato, sono ottenibili con l’istruzione:

1 F.den{1}

Per ulteriori informazioni si rimanda ovviamente alla guida di Ma-tlab.

freqrespFr

All’interno della funzione bodeFr viene richiamata la funzione ausilia-ria freqrespFr, che esegue il calcolo vero e proprio della risposta infrequenza del sistema. La sintassi è:

1 [gain, phase, w] = bodeFr(F, Q, w, delay)

in cui i parametri sono esattamente gli stessi del comando bodeFr.La funzione estrae i coefficienti dal sistema intero F, che moltiplicano

le potenze decrescenti di s. Alla variabile s viene sostituito il vettoredelle frequenze w, infine si eleva tutto il sistema all’esponente non inte-ro comune, Q. I sistema che si ottiene in questo modo, come specificatonel § 2.2.3 nella pagina 25, un sistema identificato dai multipli nα di unvalore comune di base, in questo caso al posto di α si ha la variabile Q.

nyquistFr, nicholsFr

Per calcolare il diagramma di Nyquist è disponibile l’omonima funzio-ne frazionale, la cui sintassi è

1 [gain, phase, w] = nyquistFr(F, Q, w, delay)

mentre per calcolare il diagramma di Nichols la sintassi è

1 [gain, phase, w] = nicholsFr(F, Q, w, delay)

similmente alla altre funzioni viste. Ancora una volta i parametri sonogli stessi, il funzionamento si basa sulla funzione freqrespFr

nid

Una funzione molto utile per i sistemi digitali è la funzione nid. Questapermette di trovare un’approssimazione, secondo vari metodi, dellavariabile di Laplace k sα, la sintassi è:

1 C = nid(k, v, bandwidth, n, formula, expansion, decomposition)

I primi due parametri rappresentano i termini da approssimare, cioèk * s^v. Seguono alcuni parametri che specificano l’approssimazione:bandwidth è il range di frequenze in cui deve essere valida l’approssi-mazione se il sistema è analogico, altrimenti è il periodo di campiona-mento in secondi per sistemi discreti; n è il numero di punti o l’ordinedell’approssimazione, a seconda del metodo.

La variabile formula deve essere di tipo stringa, e può assumere unodei seguenti valori:

• crone. Prima generazione dei controllori CRONE;


• carlson. Approssimazione di Carlson con almeno n zeri ed n poli;

• matsuda. Approssimazione di Matsuda con almeno n zeri e poli;

• cfehigh. Approssimazione CFE di un termine (1 + s)ν con più din poli ed n zeri, valida per frequenze maggiori di bandwidth;

• cfelow. Analoga alla cfhigh, è l’approssimazione CFE di un termi-ne(1 + 1

s)−v per frequenze fino a bandwidth;

per le approssimazioni tempo continue, oppure si deve indicare unodei seguenti valori:

• tustin. Approssimazione di Tustin;

• simpson. Approssimazione di Simpson;

• 1ofd, 2ofd, 3ofd. Approssimazioni alle differenze finite all’indie-tro rispettivamente del primo, secondo o terzo ordine;

• delta. Trasformata delta;

• impulse. Approssimazione della risposta impulsiva;

• step. Approssimazione della risposta al gradino;

per le approssimazioni discrete, tutte caratterizzate da almeno n ele-menti di ritardo.

Anche la variabile expansion deve essere una stringa predefinita, mapuò essere omessa, in questo caso sarà usato il primo metodo. Ipossibili sono:

• mcltime, mcltimeINV. Il metodo mcltime è usato per default, ilmetodo mcltimeINV è il suo inverso;

• cfe, lscfeINV. Sono i metodi per lo sviluppo in frazioni continue;

Per la teoria dei metodi di approssimazione si rimanda alle rispettivesezioni, sia per i metodi analogici § 2.4 nella pagina 27 che per i metoditempo discreti § 3.5 nella pagina 54.

L’ultima variabile in ingresso, decomposition, anch’essa opzionale,permette di scegliere se approssimare la parte intera dell’esponentev separatamente dalla parte non intera (metodo seguito per default eraccomandato) se il valore della stringa è frac, altrimenti se il valoredella variabile è la stringa all, l’ordine v è approssimato direttamente,sebbene questa procedura sia sconsigliata.

Ovviamente un comando così complesso è basato su subroutine,ciascuna delle quali può essere eseguita come script indipendente,passando gli stessi parametri (tranne formula) ad una funzione chegeneralmente ha per nome la stringa memorizzata in formula.

contfrac

Per finire, un’ultimo comando che si cita è la funzione contfrac, cherealizza l’espansione in frazioni continue di una funzione razionale,col metodo mostrato nel § 2.4.1 nella pagina 28. La sintassi è:

1 [a, b] = contfrac (x, n)

in cui x è il termine da espandere fino al grado n. I termini restituiti,a, b, identificano rispettivamente i denominatori parziali ed i numeratoriparziali della sviluppo CFE calcolato.

B I B L I O G R A F I A

[1] Al-Alaoui, M. A. (1993), «Novel digital integrator and differen-tiator», IEE Electronics Letter, vol. 29 (4), p. 376–378. (Citato apagina 62.)

[2] Al-Alaoui, M. A. (2006), «Al-Alaoui operator and the α-approximation for discretization of analog systems», Facta Uni-versitatis (NIŠ), Ser: Elec. Energ., vol. 19 (1), p. 143–146. (Citato allepagine 55 e 62.)

[3] Anatolii Aleksandrovich Kilbas, H. M. S. e Trujillo, J. J.(2006), Theory and applications of fractional differential equations,Elsevier. (Citato a pagina iii.)

[4] Balsi, M. (2001), teoria dei circuiti elettronici, Edizioni scientificheSiderea, roma. (Citato alle pagine 27 e 36.)

[5] Barbosa, R. S., Machado, J. A. e Silva, M. F. (2006), «Discre-tization of complex-order differintegrals», Second IFAC Workshopon Fractional Differentiation and its Applications, Porto, Portogallo.(Citato alle pagine 27, 28 e 59.)

[6] Barbosa, R. S. e Machado, J. A. T. (2006), «Implementationof Discrete-Time Fractional-Order Controllers based on LS Ap-proximations», Acta Polytechnica Hungarica, vol. 3 (4). (Citato apagina 62.)

[7] Barbosa, R. S., Machado, J. A. T. e Ferreira, I. M. (2005), «Pole-zero approximations of digital fractional-order integrators anddifferentiators using signal modeling techniques», Proceedings of16st IFAC World Congress. (Citato alle pagine 28, 36, 58, 59, 60

e 62.)

[8] Bohannan, G. W. (2002), «Analog realization of a fractional con-trol element – Revisited», in «Tutorial Workshop on fractional cal-culus and its applications in automatic controls and robotics»,41

st IEEE Conference on decision and control, Las Vegas, NE,USA, http://mechatronics.ece.usu.edu/foc/cdc02tw/. (Citatoa pagina 65.)

[9] Carlson, G. e Halijak, C. (1964), «Approximation of FractionalCapacitors (1/s)1/n by a Regular Newton Process», Circuit Theory,IEEE Transactions on, vol. 11 (2), p. 210–213. (Citato alle pagine 27,28 e 32.)

[10] Cervera, J., Banos, A., Monje, C. e Vinagre, B. (2006), «Tuningof Fractional PID Controllers by Using QFT», IEEE Industrial Elec-tronics, IECON 2006 - 32nd Annual Conference on, p. 5402–5407.(Citato a pagina 37.)

[11] Charef, A., Sun, H., Tsao, Y. e Onaral, B. (1992), «Fractal systemas represented by singularity function», Automatic Control, IEEETransactions on, vol. 37 (9), p. 1465–1470. (Citato alle pagine 28

e 33.)

89

http://mechatronics.ece.usu.edu/foc/cdc02tw/

90 Bibliografia

[12] Chen, Y. (2006), «Ubiquitous fractional order controls?», in «Proc.of The Second IFAC Symposium on Fractional Derivatives andApplications (IFAC FDA06, Plenary Paper.)», p. 19–21. (Citato apagina iii.)

[13] Chen, Y. Q. e Moore, K. L. (2002), «Discretization schemes forFractional-order differentiators and integrators», Circuits and sy-stems, IEEE Transactions on, vol. 49 (3), p. 363 – 367. (Citato allepagine 27, 60, 61 e 62.)

[14] Chen, Y. Q. e Vinagre, B. M. (2003), «A new IIR-type digitalfractional order differentiator», Signal Processing, vol. 83 (11), p.2359–2365. (Citato a pagina 62.)

[15] Chen, Y. Q., Vinagre, B. M. e Podlubny, I. (2003), «A new discre-tization method for fractional order differentiators via continuedfraction expansion», in «Proceedings of The First Symposium onFractional Derivatives and Their Applications at The 19th BiennialConference on Mechanical Vibration and Noise, the ASME Inter-national Design Engineering Technical Conferences & Computersand Information in Engineering Conference (ASME DETC2003)»,p. 1–8. (Citato a pagina 62.)

[16] Das, S. (2008), Functional Fractional Calculus for System Identificationand Controls, Springer. (Citato alle pagine iv, v, 40, 53, 67 e 82.)

[17] Jones, H. e Shenoi, B. (1970), «Maximally flat lumped-elementapproximation to a fractional operator immittance function», Cir-cuit Theory, IEEE Transactions on, vol. 17 (1), p. 125–128. (Citatoalle pagine 27 e 32.)

[18] Kabal, P. e Ramachandran, R. P. (1989), «Joint optimizationof linears predictors in speech coders», IEEE transaction on acou-stic speech and signal processing, vol. 37 (5), p. 642–650. (Citato apagina 80.)

[19] Kabal, P. e Ramachandran, R. P. (1989), «Pitch prediction infilter speech coding», IEEE transaction on acoustics speech and signalprocessing, vol. 37 (4), p. 467–478. (Citato a pagina 80.)

[20] Kaczorek, T. (2008), «Fractional positive continuous-Time linearsystems ant their reachability», International Journal of appliedmathematics and science, vol. 18 (2), p. 223–228. (Citato a pagina 36.)

[21] Khovanskii, A. N. (1963), The application of continued fractions andtheir generalizations to problems in approximation theory, Library ofapplyed analysis and computational mathematics, P. NoordhoffN. V., Groningen, Olanda. (Citato alle pagine 28 e 30.)

[22] Koh, C. G. e Kelly, J. M. (1990), «Application of fractional de-rivatives to seismic analysis of base-isolated models», EarthquakeEngineering & Structural Dynamics, vol. 19 (2), p. 229–241. (Citatoa pagina 78.)

[23] Laakso, T. I., Välimäki, V., Karjalainen, M. e Laine, U. K.(1996), «Splitting the unit delay [FIR/all pass filters design]», Si-gnal Processing Magazine, IEEE, vol. 13 (1), p. 30–60. (Citato allepagine 39, 43, 56, 80, 81 e 82.)

Bibliografia 91

[24] Luise, M. e Vitetta, G. M. (2003), Teoria dei segnali, McGrow-Hill,seconda edizione ed. (Citato alle pagine 13 e 43.)

[25] Machado, J. A. T. (2001), «Discrete-time fractional-order control-lers», Fractional Calculus and Applied Analysis, vol. 4 (1), p. 47–66.(Citato a pagina 62.)

[26] Monje, C. A., Calderon, A., Vinagre, B. e Feliu, V. (2004),«The fractional order lead compensator», Computational Cyberne-tics, 2004. ICCC 2004. Second IEEE International Conference on, p.347–352. (Citato a pagina 37.)

[27] Oldham, K. B. e Spanier, J. (1974), Theory and applications of diffe-rentiation and integration to arbitrary order, The fractional calculus,Academic Press. (Citato alle pagine 1, 3, 11, 12, 51 e 68.)

[28] Ortigueira, M. D. (2000), «Introduction to fractional linear sy-stems. Part 1: continuous-time case», IEE procedings - Vision imageand signal processing, vol. 147, p. 62–70. (Citato alle pagine 15, 16,25, 26 e 43.)

[29] Ortigueira, M. D. (2000), «Introduction to fractional linear sy-stems. Part 2: discrete-time case», IEE procedings - Vision image andsignal processing, vol. 147, p. 71–78. (Citato alle pagine 39, 42, 43,44, 46, 47 e 48.)

[30] Ortigueira, M. D. (2008), «An introduction to the fractionalcontinuous-time linear systems», IEEE circuit and systems magazine,p. 19–26. (Citato alle pagine iii, 7, 18, 25, 43, 54 e 82.)

[31] Oustaloup, A., Levron, F., Mathieu, B. e Nanot, F. M. (2000),«Frequency-band complex noninteger differentiator: characteriza-tion and synthesis», Circuits and Systems I: Fundamental Theory andApplications, IEEE Transactions on, vol. 47 (1), p. 25–39. (Citato allepagine 27 e 35.)

[32] Oustaloup, A., Cois, O., Lanusse, P., Melchior, P., Moreau,X. e Sabatier, L. (2006), «The CRONE approach: theoreticaldevelopments and major applications», in «Proocedings del 2°IFAC Workshop on fractional Differentiation and its applications»,Porto, Portogallo. (Citato alle pagine ix e 83.)

[33] Oustaloup, A., Cois, O., Lanusse, P., Melchior, P., Moreau, X.e Sabatier, J. (2006), «The CRONE approach: theoretical develop-ments and major applications», in «Proceedings of the 2nd IFACWorkshop on Fractional Differentiation and its Applications», .(Citato a pagina 38.)

[34] Podlubny, I. (1999), «The Laplace transform method for li-near differential equations of the fractional order», San Diego,p. 368, URL http://www.citebase.org/abstract?id=oai:arXiv.org:funct-an/9710005. (Citato alle pagine 18, 19, 25 e 56.)

[35] Roy, S. C. D. (1967), «On the realization of a constant-argumentimmittance or fractional operator», Circuit Theory, IEEE Tran-sactions on, vol. 14 (3), p. 264 – 274. (Citato alle pagine 27

e 67.)

http://www.citebase.org/abstract?id=oai:arXiv.org:funct-an/9710005

http://www.citebase.org/abstract?id=oai:arXiv.org:funct-an/9710005

92 Bibliografia

[36] Samko, S. G., Kilbas, A. A. e Marichev, O. I. (1993), FractionalIntegrals and Derivatives: Theory and Applications, Gordon & BreachScience Publishers. (Citato a pagina 1.)

[37] Schmidt, V. H., Drumheller, J. E. e Howell, F. L. (1971), «Die-lectric properties of lithium hydrazinium sulfate», Phys. Rev. B,vol. 4 (12), p. 4582–4597. (Citato a pagina 66.)

[38] Smith, J. O. (2009), Introduction to Digital Filters with Audio Applica-tions, http://ccrma-www.stanford.edu/~jos/filters/Formant_

Filtering_Example.html, online book. (Citato a pagina 80.)

[39] Spiegel, M. R. (1993), Variabili complesse, Etas libri. (Citato allepagine 15, 41 e 45.)

[40] Uncini, A. (2008), «Filtraggio adattativo e circuiti intelligenti»,Appunti del corso di circuiti e algoritmi per l’elaborazione deisegnali. (Citato a pagina 46.)

[41] Valério, D. (2005), «Ninteger v. 2.3 Fractional control tool-box for Matlab», http://web.ist.utl.pt/duarte.valerio/ninteger/Manual.pdf. (Citato alle pagine 83 e 84.)

[42] Valério, D. e da Costa, J. S. (2004), «Ninteger: a non-integercontrol toolbox for Matlab», Fractional derivatives and applications.(Citato alle pagine 24, 83 e 84.)

[43] Vinagre, B. M., Podlubny, I., Hernández, A. e Feliu, V. (2000),«Some approximations of fractional order operators used in con-trol theory and applications», Fractional Calculus and Applied Ana-lysis, vol. 3 (3), p. 231–248, http://mechatronics.ece.usu.edu/foc/cdc02tw/. (Citato alle pagine 27, 28 e 31.)

[44] Wall, H. S. (1948), Analytic theory of continued fractions, The uni-versity series in higher mathematics, D. Van Nostrand company,new York. (Citato alle pagine 28, 29, 32, 69 e 70.)

[45] Weisstein, E. W. (2009), «Binomial Coefficient», http://mathworld.wolfram.com/BinomialCoefficient.html, fromMathWorld–A Wolfram Web Resource. (Citato alle pagine 2, 3

e 58.)

[46] Weisstein, E. W. (2009), «Gamma Function», http://mathworld.wolfram.com/GammaFunction.html, from MathWorld–A WolframWeb Resource. (Citato a pagina 1.)

[47] Weisstein, E. W. (2009), «Mittag-Leffler Function», http://mathworld.wolfram.com/Mittag-LefflerFunction.html, fromMathWorld–A Wolfram Web Resource. (Citato alle pagine 18

e 19.)

[48] Xue, D. e Chen, Y. Q. (2002), «A comparative introduction of fourfractional order controllers», Intelligent Control and Automation,2002. Proceedings of the 4th World Congress on, vol. 4, p. 3228–3235,IEEE. (Citato alle pagine 36 e 37.)

http://ccrma-www.stanford.edu/~jos/filters/Formant_Filtering_Example.html

http://ccrma-www.stanford.edu/~jos/filters/Formant_Filtering_Example.html

http://web.ist.utl.pt/duarte.valerio/ninteger/Manual.pdf

http://web.ist.utl.pt/duarte.valerio/ninteger/Manual.pdf



http://mathworld.wolfram.com/BinomialCoefficient.html

http://mathworld.wolfram.com/BinomialCoefficient.html

http://mathworld.wolfram.com/GammaFunction.html

http://mathworld.wolfram.com/GammaFunction.html

http://mathworld.wolfram.com/Mittag-LefflerFunction.html

http://mathworld.wolfram.com/Mittag-LefflerFunction.html

I N D I C E A N A L I T I C O

Aanti-causale, 7, 45

armoniche, 80

BBode

diagramma, 21

Bode, 21

branch cut, 41

branch point, 41

Ccampionamento variabile, 40

capacità frazionaria, 24, 65

capacità frazionaria, viiCauchy

formule integrali, 45

formula integrazione mul-tipla, 7

Cauchy, 7

causale, 45

causalità, 7

cell arrays, 85

classe, 84

coefficienti binomiali, 2, 76

condensatore frazionario, 65

controllorePID, 37

fpid, 37

lead-lag, 36

convoluzione, 13

proprietà della, 13

costante dielettrica, 66

Dderivata, 4

differenze frazionarie, 44

differintegrale, viiGrünwald-Letnikov, 4, 6

Courant-Hilbert, 7

esponenziale, 9

potenza, 8

Riemann-Liouville, 6

discretizzazioneAl-Alaoui, 55

Bilineare di Tustin, 55

Eulero G-L, 55

metodi diretti, 54

metodi indiretti, 54

DTFT, 46

Eequazione alle differenze

intere, 44

non intere, 44

errore quadratico, 56

Eulero, 1

FFaà di Bruno, formula di, 12

filtroadattativo, 80

formant filter, 80

Fourier, ixfrattanza, 23, 24, 65, 67

frattanza, viifrattore, viifrazioni continue, 28

approssimante n-esimo, 29

denominatore parziale, 29

forma di Rogers, 28

numeratore parziale, 29

quoziente parziale, 29

Stieltjes, 28

funzionea più valori, 15

composta, differintegrazio-ne, 11

di M-L, 18

esponenziale, 18

Gamma, 1

generatrice, 55

olomorfa, 41, 45

trasferimento, 21

GGauss, 31

Grünwald, 4

HHôpital, 3

Iimpedenza, viiinduttore frazionario, 65, 67

induttore frazionario, vii

93

94 Indice analitico

integraledefinizione, 5

interpolatore, 43

LLaurent, 42

least square, 56

Leibnizformula di, 11

Leibniz, 3

Letnikov, 4

linea di taglio, 16

linearità, proprietà di, 9

linee di taglio, 41

MMacLaurin, 42

Matlab, comandibode, 24

bodeFr, 85, 86

contfrac, 87

cplxroot, 17

crone, 83

df, 73

freqrespFr, 86

gamma, 3, 53

get, 85

mlf, 18

nid, 86

ninteger, 83

roots, 79

size, 84

tf, 84, 85

vg, 73

metodo, 84

Mittag-Leffler, 18

NNichols

diagramma di, 86

norma, 56

Nyquist, diagramma di, 86

Ooggetto, 84

funzione di trasferimento,85

overloading, 84

Ppitch predictor, 80

Pochammer, simboli di, 31

polinomio indiciale, 25, 79

Post, 4

Rramo (di una funzione), 16

rapporto incrementale, 4

rete anticipatrice, 36

Riemannsuperficie di, 15, 16, 41

superficie principale, 16

Riemann, 4

risposta in frequenza, 21, 45

Sscaling, 40

Serie ipergeometriche, 31

Shannon-Whitaker, interpola-zione, 43

sistemifrazionari, viiordine non intero, vii

stabilità BIBO, 26

Stieltjes, 28

TTaylor, 42

teoremafondamentale del calcolo,

9

Parseval, 57

valore finale, 15

valore iniziale, 15

tono, 81

trasformata Z modificata, 45

trasformata di Laplace, 12

WWeil, differintegrale di, 7

ZZ

antitrasformata, 41

trasformata modificata, 46

Zadeh, 41

Circuiti lineari frazionali

Documents

Transcript of Circuiti lineari frazionali