. Algoritmii Rapizi Pentru Tratarea Semnalelor

CUPRINS

1. Introducere

Dezvoltarea algoritmilor eficienți pentru rezolvarea unor anumite probleme cheie ale tratării semnalelor a fost cu siguranță crucială pentru progresul temei noastre și pentru numeroasele sale aplicații. Algoritmul central este acela al transformării rapide Fourier, redescoperit în 1965, acesta reducând complexitatea de la O la O, ceea ce reprezintă două ordine de mărime pentru N=1024.

Analiza Fourier a devenit de atunci un instrument practic. Algoritmul este imediat aplicabil problemei convoluției obținându-se astfel câștiguri substanțiale în cazul filtrajului discret. Au mai fost dezvoltați și alți algoritmi după aceea, ameliorând performanțele în cazul spectrelor-Fourier sau în cazul filtrelor discrete.

O altă problemă importantă e soluționarea sistemelor de ecuații Toeplitz care apar în problemele este predicție liniară. Și acolo, reducerea de la O la O obținută datorită algoritmului Levinson a fost crucială pentru producerea unor calcule practice.

După o scurtă trecere în revistă a principalelor transformate, a algebrei și a unor noțiuni de teorie a numerelor, vom discuta despre diverși algoritmi de convoluție rapidă. În plus, din algoritmii optimali, bazați pe produse polinomiale, derivăm metode practice, printre altele pentru calculul convoluției liniare în cazul în care semnalul de intrare este de lungime infinită.

Apoi se va discuta despre algoritmii rapizi pentru calculul transformatei Fourier. Algoritmii Cooley-Tukey și split-radix sunt detaliați, ca și cei ai lui Winograd și a factorului prim. Pentru ultimii doi algoritmi descriem două etape importante și intermediare, permutarea Good și algoritmul Rader. Trecem apoi la algoritmii destinați în mod special transformatei multidimensionale Fourier cu generalizarea algoritmului Cooley-Tukey (numit și algoritmul vectorului radix) și a transformatei polinomiale. Niște transformate trigonometrice legate de transformata Fourier (transformata în cosinus, transformata Hartley) sunt luate în considerație, având în vedere importanța lor în aplicațiile pentru tratarea imaginii, de exemplu. Transformatele în numere întregi sunt descrise pe scurt de asemenea.

Se prezintă apoi soluția rapidă a sistemelor de ecuații cu o matrice de tip Toeplitz.

Ca urmare a acestei dezvoltări teoretice a algoritmilor rapizi, au fost discutate diverse arhitecturi, pentru tratarea semnalului precum și implementarea unor algoritmi derivați pe aceste arhitecturi.

Ca aplicație, voi prezenta filtrele Gabor și un caz particular al acestui tip de filtru, și anume filtrarea Gabor recursivă.

Acest proiect nu poate acoperi toate detaliile unui subiect atât de vast și complex, dar a încercat să prezinte într-o manieră succintă uneori importanța algoritmilor rapizi de tratare a semnalelor, importanța filtrelor Gabor și aplicațiile acestora.

2. Principalele transformate

2.1. Transformata Fourier

Transformata Fourier a unui semnal numeric este dată de

(2.1.)

unde x(K) = semnal numeric

X(f) = funcție complexă a variabilei reale continue f.

Se va utiliza și notație pentru reprezentare.

Se pune problema existenței acestei transformate.

Dacă seria converge, se spune că transformata Fourier există, și semnalul este absolut semnabil.

Pentru astfel de semnale, energia definită de

(2.2.) este o mărime finită, așadar vom vorbi de semnale de energie finită.

Fiind funcție periodică, poate fi descompusă în serie Fourier. Această descompunere nu este altceva decât transformata inversă, scrisă astfel:

(2.3.)

Integrarea se va face pe intervalul simetric față de origine .

Transformata Fourier X(f) dată de relația (2.1.) ne furnizează informații cu privire la repartiția frecvențială a amplitudinii, fazei și energiei semnalului.

Spectrele amplitudinii, fazei și energie semnalului sunt date mai jos

(2.4.)

(2.5.)

(2.6.)

Aceste formule vor conduce pentru X(f) la reprezentarea

(2.7.)

Pentru un semnal x(K) real, parte reală și partea imaginară a transformatei Fourier sunt date de:

(2.8.)

(2.9.)

Este important de remarcat că pentru un semnal real, partea reală a transformatei Fourier este o funcție pară, și partea imaginară este o funcție impară.

2.2. Transformarea Fourier discretă

Generalități

Transformarea Fourier discretă (TFD prescurtat),nu este o transformare nouă. Ea este versiunea numerică a transformării Fourier. Dacă reluăm definiția TF,

(2.10.)

constatăm că există două probleme pentru evaluarea numerică a acestei expresii. Prima este legată de variabila continuă f, pe care nu o putem introduce într-un sistem de tratare numerică, și a doua este legată de suma pe o infinitate de eșantioane, pe care nu o putem efectua în practică. Soluțiile acestor două probleme sunt imediate: trebuie să înlocuim variabila continuă f cu o variabilă discretă n și să limităm durata semnalului.

Totuși, dacă principiul soluțiilor este simplu, punerea lor în practică și interpretarea rezultatelor nu sunt simple deloc. Această secțiune este destinată studierii consecințelor transpunerii ecuației (2.10.) în caz numeric.

2.3. Discretizarea frecvenței

Variabila continuă f poate fi înlocuită cu o variabilă discretă n

(2.11.)

unde este pasul folosit pe axa frecvențelor. Frecvențele discrete se numesc frecvențe armonice. Cum x(f) este perioadă 1, este suficientă folosirea substituției (2.11.) pe o singură perioadă. Dacă N este număr de pași pe perioadă, avem:

(2.12.)

Pentru f variant de la –1/2 la +1/2, variabila discretă n variază de la –N/2 la (N/2)-1.

Ținând cont de schimbarea de variabilă, transformata inversă (2.3.) este aproximată printr-o sumă de tipul:

(2.13.)

Se folosește reprezentarea rădăcinii de ordinul N a unității prin :

(2.14.)

Astfel, (2.13.) devine:

(2.15.)

Fie valoarea exactă a acestei sume. trebuie să determinăm calitatea acestei aproximări. Trebuie să remarcăm că termenul este periodic în K de perioadă n. Astfel, semnalul , valoare exactă a sumei (2.15.) este un semnal periodic de perioadă N. Cu toate că semnalul inițial x(K) nu este presupus periodic, aproximarea obținută prin discretizarea transformatei sale Fourier este periodică. Pentru a găsi relația exactă între x(K) și , înlocuim în (2.15.) relația (2.11.) modificată cu schimbarea variabilei (2.12.). Obținem:

(2.16.)

Intercalând ordinea însumării, avem:

(2.17.)

Putem arăta ușor că expresia din paranteze face 1 și este un multiplu întreg de N și 0 în toate cazurile. Astfel:

(2.18) , pentru

Această relație indică faptul că semnalul periodic este obținut prin repetare periodică de perioadă N a semnalului x(K). Acest rezultat este dual din ceea ce decurge din teorema eșantionajului. Astfel, putem decaja următoarea regulă. Transformata Fourier, atât cea directă cât și cea inversă, a unei funcții eșantionate este o funcție periodică a cărei perioadă este inversa perioadei de eșantionaj.

Relația (2.18.) indică că, dacă durata semnalului aperiodic x(K) este limitată la N, fiecare perioadă a semnalului este o replică exactă a lui x(K).

În consecință, în acest caz, nu putem extrage x(K) exact plecând de la o perioadă . Astfel, relația aproximativă (2.15.) nu devine o identitate exactă decât pentru semnale x(K) cu durată limitată. Pentru ca astfel de semnale TFD este dată de:

(1.88) cu

Transformata inversă este dată de:

(2.19.) cu

În evaluarea relațiilor (2.19.) și (2.20.), domeniile de variație a variabilelor, respectiv K și n , sunt foarte importante. Dacă relația (2.19.) este calculată pentru toate valorile lui n, obținem un rezultat periodic. Perioada din jurul originii conduce la o reprezentare obișnuită a unor spectre frecvențiale. În schimb, dacă relația (2.20.) este calculată pentru toate valorile lui n, obținem un semnal periodic care nu are nimic în comun cu un semnal cu durată limitată decât motivul de bază. Altfel spus, dacă dispunem doar de coeficienții x(K) fără a cunoaște natura semnalului de la care provin, nu putem spune că este vorba de un semnal cu durată limitată sau periodic. Acest fapt este o ambiguitate importantă TFD. Pentru a restitui semnalul cu durată finită, este necesar să cunoaștem valoarea lui .

Este clar faptul că relația (2.10.) nu există pentru un semnal periodic. Suma nu este convergentă. Într-un asemenea caz, este util să limităm suma la o singură perioadă a semnalului. Schimbarea de variabilă (2.11.) duce în acest caz la aceleași relații (2.19.) și (2.20.) pentru semnale periodice.

2.4. Transformata în Z

Definiție

Transformata în Z a unui semnal x(K) este definită prin:

(2.21.)

unde Z este o variabilă complexă și unde x este o funcție complexă a variabilei Z. Definiția (2.21.) este numită și transformată în Z bilaterală căci sumara se extinde la toți întregii K. Ca în cazul transformatei Fourier, o serie de tipul (2.21.) pune imediat o problemă de convergență. Pentru un semnal x(K) dat, ansamblul valorilor pentru care seria (2.21.) constituie regiunea de convergență. Pentru a găsi acest domeniu folosim criteriul Cauchy pe convergența seriilor de puteri. Acest criteriu afirmă că o serie de tipul:

(2.22.)

converge, dacă următoarea condiție este satisfăcută:

(2.23.)

Pentru a aplica acest criteriu, putem descompune seria (2.21.) în două serii, una pentru K mergând (luând valori) de la – la –1 și cealaltă de la 0 la . Pentru a doua serie, condiția (2.23.) conduce la:

(2.24.)

Fie limita

(2.25.)

A doua serie converge deci pentru . Cu schimbarea de variabilă , putem arăta ușor că prima serie este convergentă pentru , unde este limita.

(2.26.)

Astfel, seria (2.21.) converge in general în un inel al planului complex al lui Z dat prin:

(2.27.)

Putem extrage din ce precede un anume număr de reguli generale pentru a găsi regiunea de convergență. Regiunea de convergență a transformatei în Z a unui semnal cauzal (nul pentru ) va fi întotdeauna exteriorul unui cerc de rază . În mod similar, regiunea de convergență a unui semnal anticazual (nul pentru ) va fi întotdeauna interiorul unui cerc de rază . Pentru un semnal cu durată finită, regiunea de convergență se constituie din ansamblul planului complex Z cu excepția punctelor Z= 0 și /sau .

Ultimul punct de remarcat este că, dacă transformata în Z este evaluată pe cercul unitate, adică dacă, transformata în Z se identifică cu transformata Fourier.

2.5. Transformata în Z inversă

Inversarea transformatei în Z poate fi stabilită aplicând teorema lui Cauchy pe integrarea pe lungul unui contor într-un plan complex. Putem decide din această teoremă că integrarea este definită prin:

(2.28.)

unde este un contur închis care înconjoară originea planului z. Multiplicând cei doi membri din relația (2.21.) prin (cu) și integrând lungul conturului care înconjoară originea și conținut în regiunea de convergență obținem:

(2.29.)

Cum integrala este calculată în domeniul de convergență, seria este convergentă, și putem interverti (inversa) integrarea și somația:

(2.30.)

Ținând cont de teorema (2.28.), obținem transformata în z inversă:

(2.31.)

Funcția de transfer a transformatei z

Să luăm mai întâi cazul unui sistem cazual. Răspunsul său impulsional este zero pentru valorile negative ale lui K. Funcția de transfer este deci dată de:

(2.32.)

Regiunea de convergență este deci interiorul unui cerc de rază . Cum o transformată în z nu converge niciodată la un pol, toți polii lui G(z) trebuie să se găsească în interiorul acestui cerc.

Dacă sistemul este stabil, dar nu neapărat cazual, atunci răspunsul impulsional este absolut sumabil:

(2.33.)

Funcția de transfer a unui sistem stabil converge în general într-un inel al planului z- urilor delimitat de .și . Ori, cum condiția (2.33.) este satisfăcută pentru un sistem stabil, funcția de transfer este convergentă deci pentru , adică pe cercul unitate. Deci, inelul de convergență trebuie neapărat să conțină cercul unitate.

3. Elemente de algebră și de teoria numerelor

3.1. Grupuri, inele, corpuri.

Să amintim pe scurt noțiunile de grupuri amintim pe scurt noțiunile de grupuri, inele și corpuri.

E vorba, în toate cele trei cazuri, de un ansamblu de elemente asociate unei (sau la două) operații și care satisface niște axiome date.

În cazul grupului, ansamblul de elemente G, asociat operației – adunare modulo 2 – (care este operația grupului și nu neapărat o adunare obișnuită), satisface:

apartenență: oricare ar fi , atunci ;

asociativitate: oricare ar fi , atunci ;

existența unei identități: astfel încât oricare ar fi ;

existența inversei: oricare ar fi astfel încât .

Precizăm că atunci când „” este comutativă, grupul este denumit abelian.

Este ușor de verificat că în orice grup, elementul identitate „e” este unic ca și inversa oricărui element al grupului. Numărul de elemente ale unui grup se numește ordinul grupului. Dacă acest număr e finit atunci e vorba despre un grup finit. În acest caz, un element compus cu el însuși și aceasta în mod repetat, va genera elementul identitate (în lumina proprietăților (i) și (iii) de mai sus). Numărul de compunere necesar pentru a genera elementul identitate pornind de la un element a se numește ordinul r al lui a.

Putem verifica că elementele succesive astfel generate formează un subgrup, și cu ajutorul descompunerii, putem arăta că ordinul oricărui element divide ordinul grupului. dacă un element g al grupului G are același ordin ca și grupul, acest element este numit generator al lui G (de vreme ce permite producerea tuturor elementelor lui G prin simpla compunere cu el însuși) grupul este calificat drept ciclic(descompunerea repetată a lui g generează cicluri).

Să trecem acum la noțiunea de inele. Adăugăm o a doua operație , ansamblului de elemente, A, și primei operații . Axiomele ce trebuie satisfăcute pentru a obține o structură a inelului sunt deci:

i) este un grup abelian;

ii) apartenență: oricare ar fi ;

iii) asociativitate: ;

distributivitate:

Dacă, în plus, operația „” este comutativă, vom vorbi de un inel comutativ. Identitatea în raport cu operația este numită 0 (zero) și dacă operația posedă un element identitate, el este unic, notat cu 1, și structura se numește inel cu identitate. Suma repetată a identității formează întregii inelului. Numărul lor e numit caracteristica inelului, care de vreme ce e finită, formează un subgrup aditiv și multiplicativ. Elementele care au inverse într-un inel cu identitate sunt numite unitățile inelului. Ele formează un grup multiplicativ. Un exemplu familiar de inel comutativ cu identitate, este mulțimea numerelor întregi cu adunarea și multiplicarea, mulțimea unităților fiind .

Mulțimea matricilor pătrate de mărime (cu elemente ) asociate adunării și multiplicării matriciale formează un inel noncomutativ, elementul identitate pentru multiplicare este matricea identitate de mărime n, iar unitățile sunt matrici nesingulare. Un alt exemplu este inelul numerelor întregi cu adunarea și multiplicarea modulo q, care este notată cu Z/(q).

O structură mai puternică decât inelul este corpul. Să luăm o mulțime de elemente C și două operații și . Structura corp este definită prin axiomele următoare:

grup abelian;

este grup abelian (0 este identitatea aditivă);

distributivitate: , oricare ar fi .

Numim 1 identitatea operației . Există deci cel puțin două elemente în orice corp: . Corpurile de mărime infinită, cum ar fi , cu operațiile de adunare și multiplicare sunt bine cunoscute. Mai puțin familiare, dar foarte interesante, sunt corpurile de mărime finite sau corpurile Galois care se scriu GF(q), unde q este număr de elemente. De exemplu, cel mai mic dintre aceste corpuri Galois este acela format de identitățile operațiilor și , adică 0 și 1 și unde și sunt adunarea și multiplicarea modulo 2. Structura algebrică a corpurilor Galois este bogată, dar explorarea sa depășește cadrul acestui capitol.

Ne vom mulțumi să menționăm că GF(q) există numai dacă q este o putere întreagă a unui număr prim, și că în acest caz, el e unic. Mai mult, orice GF(q) are un element primitiv de ordin q-1 care este generatorul grupului ciclic al elementelor diferite de zero prin operația [2].

4. Inelul întregilor

După această recapitulare a structurilor algebrice, să revenim la cazul inelului format de numere întregi și operația obișnuită de adunare și de multiplicare. Deși împărțirea nu este în general posibilă, există împărțirea cu rest. Mai exact, pentru oricare pereche de numere întregi a și b(b > 0) există o pereche de numere întregi q și r astfel încât . Numim r restul lui a în raport cu b, sau . Două numere întregi și se numesc congruente în raport cu b și . Notăm că operatorul satisface:

(4.1.)

(4.2.)

Pentru două numere întregi, cel mai mare divizor comun și cel mai mic multiplu comun sunt definiți ca cel mai mare întreg care divide a și b și respectiv că cel mai mic întreg divizibil prin a și b. Pentru a găsi cel mai mare divizor comun al lui a și b se folosește metoda următoare cunoscută sub numele de algoritmul lui Euclid.

4.1. Algoritmul lui Euclid

Luăm și luăm algoritmul și . De aici rezultă:

(4.3.)

(4.4.)

care se termină la pasul n când dă:

Cel mai mare divizor comun .

(4.5.)

Să ilustrăm metoda prin calculul CMMDC. a lui 623 și 287.

Proba constă în a nota că , oricare ar fi r și se va termina deci prin a fi 0. Astfel înmulțind matricile (nesingulare) de la 1 la n găsim:

(4.6.)

și că (ne folosim de faptul că determinantul unei matrici elementare este –1)

(4.7.)

Deci CMMDC a lui a și b divide , după (4.6.), dar divide CMMDC a lui a și b după (5.7) rezultă . Din ecuația (4.6.) rezultă:

(4.8.)

ceea ce se reduce, (atunci când) de vreme ce a și b sunt primi relativi, la identitatea lui Bézout:

(4.9.)

Vom folosi acest rezultat pentru a demonstra că aritmetica numerelor întregi modulo N (care se denotă prin ) formează un corp dacă și numai dacă N este un număr prim. Suntem îndreptățiți să verificăm că formează un inel, este suficient deci să demonstrăm existența unor inverse pentru multiplicare. Dacă N este compus și egal cu , să presupunem că are o inversă .

Ajungem la contradicția următoare:

(4.10.)

ceea ce demonstrează că nu este un corp în acest caz. Astfel, dacă N este prim, se știe că , oricare ar fi . Deci, identitatea lui Bézout ne indică că există două numere întregi și astfel încât

(4.11.)

Găsim inversa multiplicativă a lui n în felul următor:

(4.12.)

ceea ce ne arată că inversa multiplicativă a lui n este , și asta pentru oricare ar fi prim. Am arătat deci că este egal cu GF(p).

Să considerăm acum identitatea lui Bézout (4.9.), dar de această dată modulo b. Deci,

(4.13.)

Este deci clar că următoarea congruență liniară

(4.14.)

posedă întotdeauna o soluție atunci când a și b sunt numere prime relative, de vreme ce:

(4.15.)

și se adeverește că soluția astfel dată este unică. Acest rezultat poate fi folosit pentru a demonstra teorema resturilor chineze, un rezultat mult mai vechi decât algoritmul lui Euclid, dar la fel de util ca și acesta din urmă.

4.2. Teorema resturilor chineze

Dacă N se descompune în k factori primi relativi în perechi, atunci ansamblul de k congruențe liniare

(4.16.)

are o soluție unică. Ca demonstrație a teoremei, construim soluția x care satisface congruențele liniare și care este unică

(4.17.)

Pentru că este prim relativ în raport cu alți factori ai lui N, el este și prim relativ cu , și deci (4.18.) are o soluție unică după (4.15.). Această soluție poate fi calculată cu algoritmul lui Euclid în particular. Să luăm acum restul modulo a lui x în ecuația (4.17.). Toți termenii, mai puțin termenul j se anulează modulo de vreme ce este un multiplu de. Din cauza relației (4.18.), termenul j este egal cu dovedind astfel teorema.

Să discutăm pe scurt două aplicații a teoremei resturilor chineze.

Prima este aritmetica prin resturi. În loc să folosim adunarea și îmulțirea obișnuită, utilizăm în paralel niște unități aritmetice modulo evitând astfel problemele raportului. Dacă valorile de intrare și ieșire aparțin domeniului izomorfismul descris pentru teorema resturilor chineze garantează un rezultat corect. Anumite probleme practice (absența împărțirii sau diviziunii, baze de calcul exotice) au restrâns aplicabilitatea acestei metode.

O a doua aplicație a teoremei constă în rescrierea unui vector de lungime N într-un tabel cu k dimensiuni, și acesta utilizând reprezentarea:

(4.19.)

ceea ce este posibil dacă și numai dacă factorii a lui N sunt primi doi câte doi. Vom vedea prin urmare că această reprezentare multidimensională este centrală pentru dezvoltarea algoritmilor rapizi ai transformatei Fourier.

4.3. Inelul polinoamelor

O structură similară inelului întregilor cu operațiile „+” și „” este inelul polinoamelor cu coeficienți ce aparțin unui corp C și operațiile adunare și multiplicare (înmulțire) polinomiale. Să numim acest inel C(x) unde x este nedeterminată. Un polinom f(x) de grad n-1 este caracterizat prin n coeficienți ce aparțin corpului C și în care multiplică pe x. Prin extensia cazului întregilor, vom spune că un polinom A(x) divide B(x) dacă și un polinom este calificat ca fiind ireductibil dacă nu admite decât divizori de grad 0. Prin convenție, un polinom ireductibil și cu coeficientul celei mai înalte puteri egal cu 1 se numește prim. Un pol care nu e ireductibil este deci factorizabil, și această descompunere este unică. Să notăm că această factorizare depinde de corpul pe care ea este efectuată (). Ca și în cazul întregilor putem defini CMMDC și CMMMC a două polinoame a(x) și b(x).

În general, diviziunea unui polinom prin altul nu e posibilă în inelul C(x), dar putem însă evalua diviziunea cu restul, care permite scrierea, pentru două polinoame arbitrare a(x) și b(x) ().

(4.20.) , gradul

Să notăm că această reprezentare este unică. Polinomul cât q(x) și restul r(x) sunt indicați prin următoarea notație:

(4.21.)

În plus față de relațiile (4.1.) și (4.2.) care se generalizează polinoamelor, notăm că dacă c(x) este un multiplu a lui b(x), atunci:

(4.22.)

Adesea, e posibilă alegerea lui c(x) astfel încât să simplificăm calculul restului în raport cu b(x).

Algoritmul lui Euclid din paragraful prezentat anterior se aplică calculului CMMDC a lui a(x) și b(x). unde (în care) vom alege să începem (inițializăm) algoritmul astfel încât gradul .

Rezultatul etapei n pentru care este

(4.23.)

și pornind de la algoritm, identic (similar) relația (4.8), putem totuși scrie:

(4.24.)

ceea ce duce, pentru polinoamele prime relative, la relația lui Bézout:

(4.25.)

Notația zero a unui polinom, adică a unui element , astfel încât este familiară. Să ne amintim că n elemente pot fi zerouri ale lui polinom de gradul n, și că dacă este un zero a lui .

Să notăm de asemenea că un polinom nu are neapărat un zero pe corpul coeficienților săi, și că în acest caz factorizarea în termeni liniari necesită un corp de extensie potrivită. În fine, teorema resturilor chineze este în aceeași măsură generalizată cazului polinomial. Să luăm un polinom p(x) care e produsul a oricăror polinoame pi(x) prime relative în pereche. Atunci, cele k congruențe:

(4.26.)

au o soluție unică, pe care o putem construi în mod similar relațiilor (17) și (18) ca:

(4.27.)

unde sunt aleși astfel încât

(4.28.)

Proba este identică celui din cazul întregilor, adică (4.28.) are o soluție căci CMMDC și că reducerea relației (4.27.) modulo duce la relațiile (4.26.).

4.4. Produsul Kronecker

Pentru a concluziona această scurtă introducere despre instrumentele algebrice vom prezenta produsul Kronecker și una din proprietățile sale cele mai utile.

Fiind date două matrici A și B de mărime și , atunci matricea

(4.29.)

este alcătuită din blocuri de talie , în care (unde) blocul (i,j) este format de matricea .

Matricea C este deci de mărime . De exemplu cu

(4.30.)

obținem

(4.31.)

Prin substituție, putem verifica că produsul Kronecker satisface următoarea proprietate importantă (în care toate produsele trebuie să fie definite cu mărimi compatibile)

(4.32.)

4.5. Algoritmul de convoluție rapidă

Preambul

Operația de convoluție, adică filtrajul, este cu siguranță operația cea mai comună în tratarea numerică a semnalelor. Convoluția aperiodică a unui semnal a(k) cu un filtru RIF de lungime L cu un răspuns la impuls g(k) produce la ieșire b(k):

(4.33.)

Convoluția periodică sau ciclică de lungime N este obținută atunci când indicele k este luat modulo N (admitem ). Cu toate că ecuația este definită în general pentru , vom considera cazul unde a(k) este diferit de 0 pentru k de la 0 la k-1. Ieșirea e atunci de lungime k+L-1 și relația (4.31.) poate fi scrisă ca un produs de polinoame:

(4.34.) sau

(4.35.) și

Coeficienții b(i) se calculează după (4.33.). Notarea polinoamelor este nu numai mai succintă, dar și mai adecvată (potrivită) pentru dezvoltările care vin. De exemplu, convoluția ciclică se scrie cum urmează: .

Deci, convoluția periodică poate fi evaluată cu ajutorul unei convoluții aperiodice urmată de o reducere, dar în aceeași măsură și convoluția aperiodică poate fi evaluată cu ajutorul unei convoluții circulare dacă perioada satisface: . Foarte des vom evalua una din aceste convoluții înainte de a se obține pe cealaltă și invers.

4.6. Algoritmul lui Cook-Toom

Utilizăm imediat această formulare polinomială pentru a deriva un algoritm eficace de convoluție aperiodică (deci produsul polinomial). Cu ajutorul teoremei resturilor chineze, vom trece la domeniul resturilor în raport cu anumite polinoame (ale căror singură constrângere este de a fi alese două câte două). Multiplicarea polinomială este făcută separat pentru fiecare clasă de resturi, iar teorema garantează reconstrucția rezultatului pornind de la niște rezultate parțiale. Să luăm deci un polinom P(x) de grad K+L-1 în așa fel încât reducerea modulo P(x) lasă produsul neschimbat. În particular, P(x) poate fi ales ca de exemplu:

(4.37.)

Reducerea modulo constă în a evalua și rezultatul este deci de grad zero. Produsul polinomial în domeniul resturilor modulo pi(x) constă într-o multiplicare , prin modul, adică, K+L-1 multiplicări.

Reconstrucția rezultatului final de la aceste K+L-1 valori este dată de teorema resturilor chineze (4.27.) și (4.28.). Să presupunem acum că coeficienții lui A(x) și H(x), adică valorile semnalului și ale filtrului, aparțin lui , dar că zerourile ale polinoamelor utilizate în reducere sunt alese în (de exemplu numere întregi). Este deci ușor de văzut că singurele multiplicări între două elemente ce aparțin lui sunt cele K+L-1 multiplicări de . Din punct de vedere teoretic, o multiplicare printr-un element din poate fi considerată o serie de adunări (care corespund unei multiplicări datorate unui numărător ce aparține lui ) urmată de o renormalizare (împărțire prin numitor). Cum reducerea și reconstrucția conform teoremei resturilor chineze nu implică decât niște elemente ce aparțin lui , este posibil să le ridicăm la cel mai mic numitor comun (CMMNC) și să includem împărțirea prin CMMNC în multiplicare pentru. Astfel, și cel puțin teoretic, am demonstrat că, convoluția unei secvențe de lungime K cu un filtru de lungime L necesită K+L-1 multiplicări generale și nu cum ar fi în cazul unei abordări directe. Acest algoritm cunoscut sub numele Cook-Toom, poate fi în aceeași măsură interpretat ca o interpolare Lagrange, de vreme ce polinomul B(x) de grad K+L-2 este determinat în mod unic prin K+L-1 valori. Aceste valori sunt exact punctele de interpolare și la aceste puncte

Atunci când K+L este mare, acest algoritm nu e decât de interes teoretic (argumentul că multiplicarea prin numărător aparține lui devine o serie de adunări și devine imposibil de aplicat în practică în cazul mărimii acestuia). În caz contrar, pentru K și L suficient de mici, metoda produce algoritmi foarte interesanți. Vom demonstra aceasta cu un exemplu:

Pentru a calcula produsul a două polinoame de gradul doi, alegem P(x)=x(x-1)(x+1). Reducerea lui A(x) și a lui H(x) împart cu factorii lui P(x) dă:

(4.38.)

Produsul polinomial în domeniul resturilor modulo pi(x) se evaluează cu ajutorul a trei produse:

(4.39.)

Este ușor de verificat că reconstrucția se calculează astfel:

(4.40.)

Notăm deci că factorii pot fi incluși în produsul care rezultă în și și deci, avem un algoritm care necesită trei multiplicări și cinci adunări (asta în cazul în care admitem că Gi sunt precalculați – deja calculați). Un algoritm ameliorat și care nu necesită decât trei adunări poate fi derivat calculând separat și este succint dat sub forma matricială de produsul următor:

(4.41.)

Să notăm forma caracteristică a acestui algoritm. Matricea 3×2 la dreapta transformă în domeniul resturilor modulo . Matricea diagonală (formată din transformata , deja calculată) calculează produsul în domeniul resturilor. În fine, matricea 3×3 la stânga retrece valorile în domeniul coeficienților , și este deci într-un fel transformata inversă.

Această structură astfel transformată, produce punct cu punct, transformată inversă amintește calculul convoluției prin transformata lui Fourier și deci cititorul poate verifica că se obține exact această selecție atunci când elementele sunt alese ca rădăcinile de ordin N ale unității .

4.7. Algoritmi de convoluție circulară

Vom aplica exact aceeași metodă în cazul convoluției circulare. Singura diferență este că, acum, nu avem libera alegere a lui , de vreme ce el ne este dat ca conform 36. Evident, poate fi descompus pe , ceea ce va conduce la soluția obișnuită de convoluție prin transformata lui Fourier. Întotdeauna în scopul de a minimaliza numărul de multiplicări generale, vom factoriza pe , mai degrabă decât pe , ceea ce va conduce la o transformată (la fel ca și inversa sa) care nu va cere multiplicări (cel puțin în teorie). Factorizarea lui în (pe) dă niște factori care se numesc polinoame ciclotomice și se adeverește numai pentru N mic (), coeficienții acestor factori sunt toți 0, 1 sau –1.

De exemplu: se factorizează în astfel:

(4.42.)

Un algoritm de convoluție circulară de lungime 4 poate fi dezvoltat reducând și în raport cu , și . Apoi produsul resturilor este calculat întotdeauna modulo , și . Acest ultim produs rezultă din multiplicarea a două polinoame de gradul întâi, urmată de o reducere modulo . De vreme ce este primul în , folosim algoritmul lui Cook-Toom (3 multiplicări) urmat de o reducere. Am utilizat deci, toate 5 multiplicări (1, 1 și 3 pentru produsele modulo , și respectiv ) și reconstrucția este făcută din nou după teorema resturilor chineze, care nu necesită deci nici o multiplicare.

Numărul de multiplicări necesare în acești algoritmi este dat de relația următoare (funcția indică numărul de multiplicări utilizate într-un algoritm):

(4.43.)

unde k indică numărul de factori primi de . Reducerea la k multiplicări vine de la faptul că fiecare factor prim conduce la o convoluție aperiodică care necesită multiplicări, fiind gradul factorului. Forma generală a algoritmilor de convoluție dezvoltați mai sus este arătată în figura 4.1. Dacă polinoamele sunt de grad mai mare ca 1, produsul modulo este evaluat cu algoritmul lui Cook.Toom. Rezultă că algoritmii de convoluție pot fi puși, la fel cum am văzut în relația 4.41., sub forma:

(4.44.)

unde x și y sunt vectorii de intrare și ieșire, A și C matricile de adunări care corespund operațiilor de reducere și de reconstrucție, iar B este o matrice diagonală care conține reducerea (precalculată) a filtrului g. Dimensiunea (mărimea) lui B este egală cu numărul de multiplicări. Algoritmii de convoluție dezvoltați după metodele indicate mai sus și puși sub forma (4.44.), sunt des calificați ca fiind algoritmi de convoluție Winograd. Un rezultat interesant este că numerele de multiplicări obținute în algoritmul lui Cook-Toom pentru convoluția aperiodică și algoritmul lui Wingrad pentru convoluția circulară de lungime N (4.43.) sunt de fapt numerele minime posibile pentru aceste probleme. Proba acestor borne inferioare se bazează pe dimensiunea diferenței dintre spațiul vectorial al ieșirilor al unui calcul dat cu spațiul vectorial al datelor și acolo scalarii sunt restrânși la (corpul constantelor). În acest sens, algoritmii pe care îi vom dezvolta sunt niște algoritmi optimali.

Fig.4.1. Produs de 2 polinoame și modulo, un polinom având factorii . Domeniul problemă a lui este de obicei precalculat. Să notăm că modulo și că ( sunt primi relativi).

4.8. Algoritmi de convoluție lineară

Cititorul va găsi poate explicațiile de mai jos puțin prea teoretice, dar vom demonstra că algoritmii astfel dezvoltați au aplicații imediate. Mai târziu, vom utiliza aceeași algoritmi pentru a dezvolta niște metode pentru transformata Fourier foarte performante.

Am văzut că metoda lui Cook-Toom își pierdea interesul în momentul în care polinoamele nu mai erau de grad redus. O idee foarte simplă permite totuși obținerea unor câștiguri considerabile pentru produsele polinomiale de talie mijlocie. Să admitem că și au același grad impar, N-1 și să rescriem aceste polinoame ca:

(4.45.)

(4.46.)

unde conține coeficienți de puteri pare, și cei de puteri impare (și la fel pentru ). Produsul este echivalent cu produsul polinomial pe care l-am calculat pentru algoritmul descris în ecuația (4.41.), doar că coeficienții sunt acum ei însuși polinoame. Folosind (4.41.), trei produse de polinoame de grad sunt necesare, ceea ce este o reducere de 25% în numărul de multiplicări.

Dar metoda devine de-a dreptul interesantă dacă ea poate fi repetată de exemplu, dacă N este o putere de 2, . Deci, numărul de multiplicări e dat de :

(4.47.)

deci, o reducere substanțială față de multiplicări din algoritmul direct. Numărul de operații pentru produsul polinoamelor până la N=64 e dat în tabelul 4.1. Aceste metode sunt bune (dau randament) până la mărimi mijlocii, moment în care metodele prin transformată au întâietate.

Tabelul 4.1

Calculul produselor polinomiale în funcție de algoritmii mici

Foarte des, semnalul de intrare într-un filtru este de lungime infinită. Pentru a obține niște reduceri al cumulului de calcul, începem prin a evalua niște eșantioane de ieșire pe rând. De exemplu, dacă vom calcula două eșantioane dintr-un filtru de lungime 2, avem:

(4.48.)

Vom remarca că matricea de mărime este pur și simplu transpusa celei care apare într-un produs polinomial . Factorizarea eficace pentru această problemă dată în (4.41.) va da, o dată transpusă, un algoritm eficace pentru (4.48.). Din cauza formei speciale a algoritmilor de convoluție în (4.44.) care concentrează toate multiplicările într-o matrice diagonală, numărul acesteia rămâne neschimbat prin transpunere. Putem calcula deci două ieșiri din filtru în (4.48.) cu ajutorul a trei multiplicări, și aceasta în felul următor:

(4.49.)

Să notăm că numărul de adunări nu rămâne neschimbat în general în momentul transpunerii. Este interesant de arătat că acest algoritm poate fi aplicat în cazul semnalelor infinite și figura 4.2. indică schematic cum convoluția printr-un filtru de lungime 2L poate fi înlocuită și aceasta doar cu ajutorul adunărilor, prin trei filtre de lungime L și care funcționează cu viteză înjumătățită.

Aceasta concluzionează pentru moment discuția noastră despre algoritmii de convoluție rapidă. Am omis însă un anumit număr de rezultate, printre care și algoritmul lui Agrawal-Cooley care transformă o convoluție ciclică de perioadă , și primi relativi într-o convoluție bi-dimensională , cu ajutorul unui reindexaj bazat pe teorema resturilor chineze. Această metodă este din punct de vedere conceptual similară permutării lui Good pe care o vom discuta în detaliu pentru transformata lui Fourier.

În final, amintim că am dezvoltat mai ales niște algoritmi de lungime mică, care se dovedesc utili pentru algoritmii transformatei lui Fourier în special.

Fig.4.2. Calculul convoluției unui semnal de intrare infinit cu un filtru de lungime 2L cu ajutorul a trei filtre de lungime L și care funcționează cu viteza înjumătățită. Simbolurile și supraeșantionare și subeșantionare printr-un factor 2. Transformata în z a filtrului se scrie astfel:

Pentru convoluții de lungime mare, metoda clasică bazată pe proprietatea de convoluție a transformatei lui Fourier și care necesită calculul a două transformate și a N multiplicări în domeniul frecvențial (pentru o convoluție circulară de lungime N) rămâne metoda de alegere. Dar pentru a demonstra aceasta, trebuie mai întâi să dezvoltăm niște algoritmi rapizi pentru transformata lui Fourier.

4.9. Algoritmi ai transformatei Fourier.

Generalități

Ideea de la baza tuturor algoritmilor transformatei rapide Fourier (TFR) e aceea de a împărți problema inițială în sub-probleme și aceasta cu un cost cât mai scăzut posibil. Această metodă, cunoscută sub numele „împarte și cucerește”, poate fi repetată în general pe aceste sub-probleme și aceasta conduce la niște reduceri impresionante a numărului de calcule:

(4.50.)

Notăm că algoritmul final este în general foarte sensibil la costul transformării și la costul sub-problemelor minime (acelea care nu mai pot fi subdivizate). Trei clase principale de TFR vor fi considerate, și aceasta în funcție de lungimea N a transformatei:

N compus, cu factori arbitrari

N compus, cu factori coprimi

N prim.

Algoritmul lui Cooley-Tukey

Vom începe de la cazul cel mai cunoscut și redescoperit de Cooley-Tukey și care a fost la baza expansiunii explozive a aplicațiilor tratării numerice a semnalelor. Amintim că transformata Fourier (TF) este definită de perechea:

(4.51.)

(4.52.)

unde este rădăcina de ordin N a unității ().

Algoritmul Cooley-Tukey se aplică de fiecare dată când N este compus, dar dacă factorii lui N sunt primi relativi, vom descrie o metodă mai performantă datorată lui Good. Rezultă ca algoritmul Cooley-Tukey este des utilizat pentru lungimi care sunt puteri ale unui număr prim (și în general mic, cum ar fi 2 sau 3).

Înainte de a deriva forma generală a algoritmului, să considerăm cazul și să derivăm un algoritm în baza 2, care consistă în a împărți mărimea problemelor prin 2 la fiecare etapă. Să luăm de exemplu n pereche, în (4.51.).

(4.53.)

Am folosit faptul că . În (4.53.) avem o rădăcină de ordin N/2 a unității (care corespunde unei transformate Fourier de lungime N/2) dar suma se efectuează pe N termeni. Din contra, și sunt multiplicați de același factor căci:

(4.54.)

deci, (4.53.) se rescrie după cum urmează:

(4.55.)

adică termenii cu index par a transformatei Fourier de lungime N sunt obținuți cu ajutorul unei transformate Fourier de lungime N/2 pe o secvență formată prin:

Cazul impar duce la:

(4.56.)

Să notăm că și sunt multiplicate de (prin) același factor în timp ce:

(4.57.)

Urmează deci că:

(4.58.)

adică că termenii cu indici impari ai transformatei Fourier de lungime N sunt obținuți pornind de la o transformată Fourier de lungime N/2 pe o secvență formată de .

Deci, cu prețul a N/2 multiplicări, am redus problema unei transformate Fourier de lungime N la cea a două transformate Fourier de lungime N/2. presupunând că o transformată Fourier de lungime N folosește multiplicări, noi verificăm că inegalitatea (4.50.) este verificată de vreme ce:

(4.59.)

Evident, dacă N este o putere de 2, metoda poate fi repetată, și după aplicații, sub-problemele sunt transformate Fourier de lungime 2, adică sume și diferențe. Numărul de multiplicări complexe deci la, o complexitate mai mică decât cele multiplicări utilizate în algoritmul direct.

Reconsiderând (4.56.) și (4.58.), se poate spune că o transformată Fourier de lungime 2 se aplică perechilor de eșantioane și , că apoi, o serie de multiplicări prin puteri succesive de se aplică diferenței și că apoi două transformate Fourier de lungime N/2 sunt evaluate pornind de la secvențe intermediare. Această structură e generală în algoritmii Cooley-Tukey, după cum vom vedea.

Dacă , putem rescrie indicii n și k astfel:

(4.60.)

(4.61.)

Această metodă nesimetrică de a reprezenta indicii va conduce la o simetrie la nivelul transformatei Fourier. Putem verifica că această reindexare este completă. Utilizând această formă pentru n și k în transformata Fourier în (4.51.) obținem:

(4.62.)

Să notăm că:

(4.63.)

Deci, (4.62.) poate fi scris astfel:

(4.64.)

de la dreapta la stânga:

N1 transformata Fourier de dimensiune N2 ,

în jur de multiplicări prin puteri ale rădăcinii de ordin N a unității (numită twiddle factors în engleză),

N2 transformata Fourier de dimensiune N1

Fig.4.3. Calculul unei transformate Fourier (TF) de lungime prin algoritmul Cooley-Tukey

Ca exemplu, figura 4.3. indică schematic calculul unui T.F. de lungime .

Un alt exemplu este acela a unei transformate Fourier în baza 4.

Dacă , putem rescrie (4.64.) sub formă matriceală după cum urmează:

(4.65.)

sau

Structura algoritmului apare deci clar, cu cele 4 transformate Fourier de lungime N/4 urmate de multiplicarea punct cu punct și de o transformată Fourier de lungime 4 în ieșire.

Notăm că factorizarea dă loc unei structuri similare dar care începe printr-o T.F. de lungime 4 și se termină printr-o T.F. de lungime N/4, aceste două forme ale algoritmului se numesc decimarea în timp și în frecvență, de vreme ce (4.65.) folosește versiuni decimate ale secvenței de intrare, pe când cealaltă formă produce versiuni decimate ale ieșirii.

Algoritmul cu bază mixtă

Să luăm acum algoritmii de decimare în frecvență cu baza 2 și 4 și asta pentru transformate ale căror lungimi sunt puteri de 2. Avantajele relative ale celor 2 algoritmi sunt următoarele: algoritmul în bata 2 calculează partea pară (5.5) a ieșirii cu un cost multiplicativ nul; algoritmul în baza 4 nu calculează decât o pătrime din ieșiri, în ocurență , cu cost multiplicativ nul, din contră, numărul de etape în algoritmul este aproape de 2 ori mai mic decât algoritmul în baza . Costul multiplicativ total fiind egal cu numărul etapei ori cu numărul de multiplicări din fiecare etapă, se vede că acești doi algoritmi reprezintă un compromis diferit între numărul și complexitatea etapelor algoritmului. Este deci normal să considerăm o abordare hibridă care să lege avantajele celor doi algoritmi ceea ce este făcut în algoritmii cu bază mixtă (sau split-radix in limba engleză). Pe scurt, parte cu indici pari a ieșirii este evaluată cu baza 2, iar restul este calculat cu baza 4. Aceasta dă naștere următoarelor ecuații:

(4.66.)

(4.67.)

(4.68.)

Deci o T.F. de lungime N a fost redusă la o T.F. de lungime N/2 și la două T.F. de lungime N/4, și asta cu prețul a N/2 multiplicări complexe prin rădăcini de ordinul N ale unității și N adunări complexe. Se dovedește că acesta este cel mai bun compromis posibil în termeni de numere de operații. Figura 4.4. arată schematic evaluarea unei T.F. de lungime 16 cu ajutorul algoritmilor cu baza 2, 4 și respectiv de bază mixtă.

În termeni de complexitatea calculului, o numărare a numerelor de operații (adunări și multiplicări) la fiecare etapă conduce la niște formule recursive care, combinate cu complexitățile problemelor inițiale (T.F. de lungime 2 și 4) dau complexitatea pentru o T.F. de lungime oarecare. Să notăm că o multiplicare complexă folosește fie sau (unde semnifică (înseamnă) multiplicare respectiv adunare reală) și că vom alege algoritmul cu 3 multiplicări. Apoi, înmulțirea prin nu utilizează operații aritmetice iar înmulțirea cu (sau cu una din puterile sale impare) folosește. În fine, o T.F. de lungime 2 folosește și cea de lungime 4 folosește . Ținând cont de aceste numere, verificarea că o T.F.R. cu baza 2 folosește (N o putere a lui 2):

(4.69.)

(4.70.)

este posibilă.

O T.F.R. cu baza 4 conduce la (N o putere pară de 2):

(4.71.)

(4.72.)

Și, în fine, T.F.R. cu bază mixtă cere (N o putere de 2):

(4.73.)

(4.74.)

Fig.4.4. Calculul unui T.F. de lungime 16 cu ajutorul unor diferiți algoritmi: (a) algoritm cu baza 2; (b) algoritm cu baza 4; (c) algoritm cu bază mixtă (split-radix).

Asimptotic, algoritmii cu baza 2 și 4 (solicită) cu 50% și 12,5% mai multe multiplicări decât algoritmii cu bază mixtă. Tabelele 4.2. și 4.3. indică numărul de multiplicări și de adunări ale acestor trei algoritmi pentru lungimi de la 16 la 2048. Evident, se pune întrebarea dacă algoritmul cu bază mixtă este optimal. Din punctul de vedere al complexității multiplicative, este posibil să obținem algoritmi cu structură similară, dar mai eficace grație utilizării produsului polinoamelor optimale. Acești algoritmi nu sunt de fiecare dată aplicabili în practică, având în vedere numărul astronomic de adunări pe care îl cer.

(4.75.)

Pentru moment, din punct de vedere practic, algoritmul în bază mixtă atinge cel mai mic număr cunoscut de operații.

Tabelul 4.2.

Numărul de multiplicări pentru algoritmul în bază 2, 4, mixtă (SRFFT), factori primi (P.F.A.) și Winograd.

Tabelul 4.3.

Numărul de adunări pentru algoritmul în baza 2, 4, mixtă (SRFFT), factori primi (PFA) și Winograd.

Algoritmul Rader

Algoritmii prezentați mai sus se aplică în cazul în care N este compus. La cealaltă extremă se situează cazul în care N este prim. La prima vedere, acest caz pare un pic mai greu de tratat, pentru că matricea D.T.F. nu pare să aibă o factorizare evidentă cum a fost cazul pentru N egal cu o putere a lui 2 de exemplu. Cu toate acestea, doi algoritmi datorați lui Bluestein și a lui Rader demonstrează că este posibilă transformarea lui T.F. într-o convoluție de lungime asemănătoare. Acest rezultat este important din punct de vedere conceptual, de vreme ce, am avea de gând mai degrabă să folosim o T.F. pentru a evalua o condiție. În plus, în cazul algoritmului lui Rader, este vorba de un algoritm foarte performant, mai ales dacă este conjugat cu niște algoritmi de convoluție optimali derivați în secțiunea 4.3.

Atunci când T.F. este de lungime N egală cu p, un număr prim, algoritmul lui Rader profită de faptul că indicii (întregii modulo p) împreună cu adunarea și înmulțirea modulo p formează un corp Galois, (GF(p)). În particular, am văzut la secțiunea 4.2. că există cel puțin un element primitiv, , încât elementul din GF(p) poate fi exprimat ca una din puterile p-1 a lui . Să numim r(n) puterea lui astfel încât

(4.76.)

Deci, r(n) joacă un rol similar logaritmului în baza pe corpul GF(p). Toți indicii, pornind de la indicele 0, pot fi exprimați ținând cont de (4.76.). Să separăm deci x(0) de restul T.F. și scriem:

(4.77.)

(4.78.)

Evident, toată complexitatea multiplicativă rezidă în suma de la dreapta în (4.78.). Să numim această sumă care duce la .

Înlocuim acum toți indicii în prin notarea echivalentă, după (4.76.):

(4.79.)

Notăm că în cazul în care r(n) ia valori de la 0 la p-2, n traversează toți întregii de la 1 la p-1 (după o evidentă permutare). Putem deci rescrie suma de (4.79.) în funcție de r(n):

(4.80.)

care este o permutare ciclică între secvențele și . Dacă ordinul uneia din secvențe este inversat, obținem în mod similar o convoluție. Cum secvențele care apar în (4.80.) sunt doar niște permutări ale secvențelor de intrare și ieșire, algoritmul lui Rader permite evaluarea unei T.F. de lungime p ca o convoluție circulară de lungime . Să luăm ca exemplu o T.F. de lungime 5:

(4.81.)

Pe noi ne interesează decât submatricea din partea dreaptă jos. În GF(5), elementul 2 este un element primitiv (). Reordonând intrările și ieșirile în așa fel încât să apară indicii care corespund puterilor succesive ale lui 2 modulo 5, adică:

(4.82.)

submatricea care ne interesează în (4.48.) devine:

(4.83.)

care este o matrice de corelare circulară. Acum, este suficient să derivăm un algoritm de convoluție circulară optimal, de lungime 4, adică un produs polinomial, modulo . Aceasta dă un algoritm în forma (4.44.). Se dovedește că, de vreme ce secvența nu este arbitrară corespunde unei urmări particulare de putere a rădăcinii de ordinul p a unității, matricea diagonală B nu conține decât niște elemente pur reale sau imaginare. Aceasta conduce deci la un număr de înmulțiri (după (4.43.) și remarca precedentă) egal cu:

(4.84.)

unde k este un număr de factori de . În exemplul T.F. de lungime 5, aceasta dă loc la 10 multiplicări, a se vedea (4.42.).

Se dovedește că algoritmul Rader poate fi generalizat la niște puteri de numere prime, și așa putem găsi algoritmi optimali (în sensul complexității multiplicative), de exemplu, pentru transformate de lungime egală cu o putere a lui 2 (4.75.).

Permutarea Good

Am văzut că metoda Cooley-Tukey poate fi aplicată din momentul în care lungimea N a T.F. este compusă. Să ne amintim că costul transformării nu e nul, de vreme ce un anumit numărul de multiplicări este inevitabil. Se dovedește că, de vreme ce factorii lui N sunt primi relativi, o transformare în subprobleme poate fi obținută care să nu antreneze nici un cost în termen de complexitate aritmetică. În locul permutărilor date prin relațiile (4.60.) și (4.61.), folosim .

(4.85.)

(4.86.)

Este indicat să verificăm că (4.85.) și (4.86.) reprezintă niște permutări și că doar și sunt primi relativi. Înlocuind aceste expresii în formula T.F. (4.51.), se știe că exponentul rădăcinii de ordinul N a unității poate fi luat modulo N. Deci:

(4.87.)

în virtutea faptului că . Deci, T.F. (4.51.) devine

(4.88.)

Aceasta este o transformată Fourier bidimensională de mărime dar cu exponenți și permutați prin respectiv . Pentru a obține ordinul obișnuit al exponenților, folosim teorema resturilor chineze (4.17.) și (4.18.) prin (cu):

(4.89.)

unde,

(4.90.) și

În acest caz, (4.88.) devine:

(4.91.)

Să notăm că rolurile lui și pot fi schimbate între ele. Toată dezvoltarea descrisă mai sus se generalizează în cazul în care N are un număr oarecare de factori primi relativi (fie folosind metoda descrisă mai sus, fie folosind teorema resturilor chineze).

Algoritmul factorilor primi

Atunci când transformata bidimensională este evaluată direct, aranjând datele într-un vector cu două dimensiuni și calculând transformatele pe liniile și coloanele acestui vector, se obține algoritmul factorilor primi.

Evaluarea (4.88.) și (4.89.) cer calculul lui , și de dimensiune . Notând cu și numărul de înmulțiri și adunări ale unui T.F. de dimensiune , se dovedește că:

(4.92.)

(4.93.)

În cazul în care N are d factori primi relativi:

(4.94.)

numerele de multiplicări și adunări sunt:

(4.95.)

(4.96.)

Aceste numere de operații sunt indicate în tabelele 4.2. și 4.3. pentru lungimi de interes practic.

Figura 4.5. indică schematic calculul unei T.F. de lungime 15 pentru algoritmul factorilor primi. O comparație cu figura 5.4 pune în evidență faptul că multiplicările intermediare prezente în algoritmul Cooley-Tukey au dispărut grație permutării Good.

Fig.4.5. Calculul unei T.F. de lungime 15 cu algoritmul factorilor primi.

Algoritmul Winograd

Am văzut că permutarea Good transformă o T.F. de dimensiune într-o T.F. de lungime prin . Să numim și matricile care corespund unei T.F. de mărime și respectiv. Matricea F corespunde unei transformate bidimensionale și poate deci fi scrisă și ajutorul produșilor Kronecker ca:

(4.97.)

Vectorul de intrare este format pornind de la un tabel bidimensional (obținut după permutarea Good) prin plasarea coloanelor unele după altele. Acum, dacă și sunt calculați cu algoritmul Rader, este vorba despre matrici de convoluție. Un algoritm optimal pentru va avea deci forma caracteristică din ecuația (4.44.) în care matricea conține toate înmulțirile (multiplicările). Să notăm că dacă este de mărime sunt deci multiplicări, dar că cel puțin una din acestea e trivială (cea care corespunde calculului unei T.F. de frecvență zero, a se vedea drept exemplu ecuația (4.81.)). Această distincție între mărimea lui și numărul de elemente diagonale care corespund unor multiplicări netriviale este necesară în evaluarea complexității multiplicative a algoritmului Winograd.

Să reluăm formularea (4.97.) dar să înlocuim acum matricile prin niște algoritmi optimali:

(4.98.)

Utilizând proprietatea produșilor Kronecker indicată de relația (4.32.) rezultă că:

(4.99.)

Toate multiplicările se găsesc regrupate în matricea diagonală centrală care este de mărime . Dacă și conțin și multiplicări triviale, atunci cer

(4.100.)

unde am folosit faptul că .

Este vorba de compararea relației (4.100.) cu rezultatul echivalent pentru algoritmul factorilor primi. Ținând cont de faptul că în algoritmii optimali ai T.F. numărul de multiplicări este întotdeauna strict mai mic decât (unde este lungimea unei T.F.), putem verifica că algoritmul Winograd folosește întotdeauna mai puține multiplicări decât cel l factorilor primi.

În ceea ce privește adunările, se dovedește că numărul lor depinde de ordinul în care aplicăm cele două T.F., căci numărul lor

(4.101.)

Fig.4.6. Calculul unei T.F. de lungime 15 pentru algoritmul Winograd

Calculul unei T.F. de lungime 15 pentru algoritmul Winograd nu este simetric în și. Va fi vorba deci de alegerea ordinului care minimizează (4.101.). Să notăm de asemenea că de vreme ce algoritmul Winograd folosește cel puțin la fel de multe adunări ca algoritmul factorilor primi, a se vedea ecuația (4.103.).

În general, numărul adunărilor este net superior și tehnicile speciale(cum ar fi split-nesting) au fost dezvoltate înainte de a-l reduce.

Figura 4.6. descrie schematic calculul unei T.F., de lungime 15, de data asta cu algoritmul lui Winograd. Să notăm că este implantarea directă a formulei (4.109.). Numerele de operații pentru lungimi până la N=2520 sunt date în tabelele 4.2 și 4.3.

Este interesant să comparăm algoritmii cu bază mixtă ai factorilor primi și Winograd (pentru lungimi comparabile, de vreme ce acești algoritmi nu se aplică acelorași lungimi). Din punct de vedere al numărului de multiplicări , algoritmul Winograd este bineînțeles cel mai eficace, urmat de cel al factorilor primi și apoi de cel cu bază mixtă. În cazul adunărilor, situația este inversă. Dacă vom compara acum cu trei algoritmi în funcție de numărul total de operații (înmulțiri și adunări), vom descoperi că ele se urmează îndeaproape. Acest lucru este arătat în figura 4.7.

Până în prezent ne-am limitat să discutăm despre complexitatea aritmetică. Ne punem întrebarea ce o fi cu această complexitate structurală de influențează atât viteza de execuție a unei implantări logice cât și costul unei realizări materiale.

Fig.4.7. Numărul total de operații pentru algoritmii cu bază mixtă, pentru algoritmii primi și pentru algoritmii Winograd.

Bineînțeles, spre deosebire de complexitatea aritmetică, care are o definiție foarte precisă, această complexitate structurală rămâne relativ vagă (neclară). Este clar totuși că permutările cerute și regularitatea acesteia joacă un rol cheie, de exemplu, pentru numărul de transferuri în memorie (care poate reprezenta o parte neglijabilă a timpului de execuție a unei T.F.R.). Pe scurt, algoritmii cel mai bine structurați (cum ar fi cei în baza 2, 4, sau în bază mixtă) vor avea un avantaj cert în raport cu algoritmii cu structură complexă, cum ar fi algoritmul Winograd. Algoritmul factorilor primi este în mod cert de preferat algoritmului Winograd, dar nu e atât de simplu cum sunt algoritmii pentru puteri ale lui 2. Avem deci o situație inversă celei complexității multiplicative.

Vorbind despre implantări, să notăm că algoritmii cu excepția celui Winograd care atinge dimensiunea unui vector cu date intermediare)pot fi implantați într-un mod numit „la loc”, adică care reutilizează întotdeauna aceleași poziții din memorie (plus un număr limitat din memoria temporară). În acest caz, rezultatele ieșirilor sunt permutate. Aceasta nu are nici o importanță dacă T.F. este folosită de exemplu pentru a calcula convoluția, de vreme ce este suficient să ținem cont de această permutare în calculele ce urmează (convoluție în domeniul T.F. inversă). Din contră, dacă rezultatele trebuie să fie în ordine, este vorba de a folosi niște metode rapide pentru a inversa această permutare care este specială. De exemplu, pentru algoritmii în baza 2, ordinul este reprezentat în baza 2 dar inversând ordinul biților. Mai mulți algoritmi rapizi au fost propuși pentru a învinge această permutare.

În fine, diverși algoritmi descriși mai sus vor avea un comportament numeric diferit. Erorile de rotunjiri și propagările lor depind de numărul de operații și de structura algoritmului.

Algoritmi ai transformatelor multidimensionale

În secțiunea precedentă, am întâlnit deja o transformată multidimensională ca rezultat al permutării Good (4.88.). Am descris deci două metode de evaluare, una fiind cea a liniilor/coloanelor și dând loc algoritmului factorilor primi, celălalt dând loc algoritmului Winograd. Aceste două metode sunt la fel de posibile atunci când transformata este o transformată bidimensională adevărată (adică nu derivă dintr-o transformată monodimensională). Amintim de asemenea că în (4.88.), transformata de mărime are dimensiuni care sunt numere prime relative. Luăm aici cazul transformatei de mărime , unde N este de exemplu o putere a lui 2. În afara celor două metode precedente, putem acum utiliza două abordări noi. Prima este o extensie imediată a algoritmului Cooley-Tukey, în timp ce, celălalt mai puternic și original, se bazează pe transformate noi numite transformate polinolmiale.

Algoritmul bazei vectoriale.

Transformata Fourier bidimensională de dimensiune este definită după cum urmează:

(4.102.)

unde W este o rădăcină de ordinul N a unității.

Pentru a fi specific, alegem N ca o putere de 2 și dezvoltăm un algoritm cu bază vectorială 2. Prin extensie al algoritmului monodimensional de decimare în timp, este clar că trebuie divizate intrările și ieșirile în patru grupe, și acestea în maniera următoare:

(4.103)

(4.104)

Fiecare subansamblu conține eșantioane. Înlocuind această descompunere în definiția (4.102.), o putem rescrie pe aceasta sub formă matriceală ca:

(4.105.)

unde T.F. reprezintă o T.F. de dimensiune și unde

Remarcăm (observăm) că matricea stângă este o T.F. de dimensiune , și, deci, am redus o T.F. de dimensiune în 4 T.F de dimensiune apoi înmulțiri (multiplicări) complexe și în fine T.F. de dimensiune (fie, de vreme ce o T.F. de dimensiune folosește opt adunări, adunări)

Avem

(4.106.)

(-2 apare pentru că o T.F. folosește numai adunări).

Algoritmul în linie coloană va folosi T.F. de lungime N adică în jur de multiplicări, de unde rezultă că algoritmul în bază vectorială ameliorează performanța cu 25% în acest exemplu simplu. Evident, metoda se generalizează la baze mai mari (4,8 etc.) și tabelul 4.4. dă numărul de multiplicări prin punct pentru diverși algoritmi.

Transformatele polinomiale

Transformatele polinomiale, introduse de Nussbaumer în 1977 pentru a calcula niște convoluții multidimensionale, sunt o generalizare a transformatei Fourier. Vom considera utilizarea lor în calculul unei TF multidimensionale și asta pentru N egal cu o putere a lui 2. Rescriem (4.102.) sub formă polinomială:

(4.107.)

(4.108.)

(4.109.)

Dacă datele sunt într-un tabel bidimensional, remarcăm că (4.107.) formează un polinom cu linia , (4.108.) calculează o transformată pe fiecare coloană (reducerea mod nu are efect de vreme ce gradul lui este mai mic decât N) și în fine, (4.109.) calculează o transformată pe linii înlocuind Z cu. Cele trei ecuații de mai sus sunt deci echivalente relației (4.102.).

Acum dacă este impar este un factor de care, el însuși este un factor de . Deci, putem reduce (4.108.) modulo unde, cu

(4.110.)

(4.111.)

Cum sunt primi relativi, permutarea este completă. Înlocuim cu în (4.110.).

(4.112.)

(4.113.)

De vreme ce în (4.113.), putem înlocui cu z în (4.112.) pentru a obține:

(4.114.)

Aceasta este o transformată polinomială, cu rădăcina z definită modulo . Această transformată nu folosește multiplicări, de vreme ce polinoamele sunt multiplicate prin o putere de z (adică coeficienții sunt deplasați) apoi o reducere modulo care nu folosește adunări, este aplicată. În fine, o structură eficace de TFR poate fi folosită pentru a calcula (4.114.) cu ajutorul a adunări. După calculul relației (4.114.) e clar că relația (4.113.) nu conține decât polinoame de grad inferior lui N/2. Aceste polinoame sunt apoi evaluate pentru în virtutea relației (4.113.), adică (4.113.) corespund evaluării a N TF impare de lungime N/2 (ele se numesc impare pentru că ele corespund părții impare a unei TF de lungime N). Ceea ce este important de notat este că TF bidimensională rezultă aici într-o serie de TF monodimensionale, în orice caz pentru partea pe care am luat-o în considerare. În fine, argumente similare pot fi utilizate pentru celelalte părți în mod egal, dar trimitem cititorul interesat la referințele bibliografice pentru o dezvoltare detaliată a algoritmului complet.

Tabelul 4.4.

Numere de multiplicări prin punct pentru diverși algoritmi ai TF bidimensionale

Acest algoritm este evident mai puțin trivial decât un algoritm în linie/coloană sau cu bază vectorială, dar complexitatea rezultantă este de asemenea mult inferioară după cum putem vedea în comparația din tabelul 4.4. Din punct de vedere practic, se dovedește că programarea algoritmului prin transformată polinomială nu este ușoară, și că un algoritm în bază vectorială 4 poate fi un compromis între complexitatea aritmetică și complexitatea logică. Notăm de asemenea că un algoritm în bază vectorială 2 este mai puțin eficace decât cel mai bun algoritm monodimensional (algoritm în bază mixtă) folosit doar în linie/coloană.

4.10. Algoritmi derivați ai transformatei Fourier

Preambul

Există un număr de variații pe tema transformatei Fourier care au importanța lor pentru tratarea semnalului. O primă variație privește semnalele care sunt transformate: acestea sunt mai degrabă reale, simetrice decât complexe. Cazul real este important încât o transformată a fost dezvoltată special pentru acest caz (transformata numită Hartley). Apoi, vom descrie o aproximare a transformatei Karhumen-Loive, numită transformata în cosinus, care este folosită foarte des în codarea imaginilor, și se dovedește că ea este legată direct de transformata Fourier. În final, vom vorbi pe scurt despre transformatele în numere întregi, care sunt o generalizare a transformatei Fourier pe corpuri finite.

Transformata Fourier pentru semnale condiționate

Foarte des în aplicații, semnalele de tratat sunt reale. Este normal să se aștepte la o dimensiune a complexității de ordinul 50%, de vreme ce jumătate din termenii de ieșire sunt redundanți din cauza simetriei ermetice a transformatei Fourier a unui semnal real.

(4.115.)

Vom deriva algoritmul cu bază mixtă pentru cazul semnalelor reale și vom demonstra astfel reducerea în complexitate. Reluând relațiile (4.66.), (4.67.), (4.68.), notăm că prima ecuație reprezintă o TF de lungime N/2 pe un semnal real:

Apoi, (4.67.) este o TF de lungime N/4 pe un semnal complex ca și N/4 multiplicări complexe. În fine, notăm că (4.68.) devine redundantă căci:

(4.116.)

Deci, și sunt complexe conjugate. Complexitatea algoritmului asfel obținut, care calculează o TF reală cu ajutorul unei TF reale de lungime înjumătățită și a unei TF complexe cu un sfert de lungime și a N/4 multiplicări complexe, este de:

(4.117.)

(4.118.)

Comparând cu complexitatea în cazul complex în (4.73.) și (4.74.), vedem că complexitatea multiplicativă a fost exact redusă la jumătate și că complexitatea aditivă este chiar un pic mai mult. O reducere similară este posibilă în cazul celorlalți algoritmi (baza 2, 4) ca și în algoritmii factorilor primi și în cel al lui Winograd. Pentru aceste două ultime cazuri, amintim că în acești algoritmi toate multiplicările sunt regrupate în matrici diagonale, și că coeficienții în aceste multiplicări sunt fie pur reali, fie pur imaginari. Atunci când intrările sunt reale, trebuie mai degrabă una decât două multiplicări pe element diagonal, și avem deci o reducere de 50% a complexității multiplicative. Notăm că structura algoritmilor reali este mai complicată decât cea a echivalenților complecși

Transformata Hartley

Această transformată este apropiată de transformata Fourier. Definiția se este suma părții reale și imaginare a unei transformate Fourier, deci:

(4.119.)

Din cauza ortogonalității sinusului și cosinusului, transformata Hartley este propria sa inversă.

(4.120.)

Evident, convoluția în domeniul temporal se simplifică în domeniul frecvențial, dar nu se reduce la multiplicări punct cu punct (având în vedere că transformata Fourier este matricea vectorilor proprii ai matricei de convoluție). Forma convoluției în domeniul Hartley este:

(4.121.)

unde și corespund părții pare și impare a transformatei Hartley a filtrului (care sunt și parte reală și imaginară a transformatei Fourier a filtrului). Relația (4.121.) se deduce ușor din proprietatea convoluției a transformatei Fourier. Notăm că în domeniul Hartley convoluția folosește 3/2 multiplicări prin punct dacă și sunt evaluați simultan, dat fiind că:

(4.122.)

Produsul matrice-vector de mai sus corespunde unei multiplicări complexe (folosind deci trei multiplicări reale), ceea ce demonstrează echivalența între această metodă și cea bazată pe transformarea Fourier. În fine, o convoluție de secvențe reale calculate cu transformate Fourier optimizate pentru acest caz folosește exact același număr de multiplicări și în mod lejer mai puține decât o metodă bazată pe transformata Hartley. Deci, singurul avantaj a transformatei Hartley este simplicitatea conceptuală a unei transformate care este propria sa inversă.

Transformata în cosinus

În aplicări cum e codarea imaginilor, locurile de eșantionare sunt codate împreună. Pentru a minimiza redundanța trebuie să decorelăm aceste eșantioane înainte de a face codajul propriu zis. Transformata optimală este cea Karkunen-Loive, ai căror vectori sunt vectorii proprii ai matricii de autocorelare a eșantioanelor. Să numim x și y vectorii înainte și după transformata T.

(4.123.)

și să alegem liniile lui T ca vectori proprii ai matricii de autocorelare R.

(4.124.)

atunci

(4.125.)

unde este matricea diagonală a propriilor valori. Am folosit deci faptul că R este simetric (deci T este ortogonală dacă normalizăm vectorii proprii. Dacă recunoaștem că eșantioanele provin dintr-un proces markovian de ordin prim (ceea ce este o aproximație grosieră dar suficientă), atunci R are următoarea formă:

(4.126.)

În plus, dacă factorul de corelare p este relativ mare (de exemplu), atunci se dovedește că transformata următoare, numită transformata în cosinus, este o aproximare bună a transformatei Karkunen-Loive.

(4.127.)

Pentru a obține o transformată unitară trebuie să multiplicăm cu . Notăm că această transformată rezultă din două aproximări succesive: o aproximare a procesului stocastic ale cărui eșantioane sunt soluții, apoi o aproximare a transformatei Karkunen-Loive pentru acest proces.

Se dovedește că transformata în cosinus, deși apare aproape dublă ca aproximare se poate aplica cu succes în aplicațiile codării, cum ar fi compresia imaginilor. Pentru partea noastră, ne vom ocupa de algoritmii eficace pentru calcului relației (4.127.).

Mai mulți algoritmi sunt posibili și vor depinde de lungimea transformatei. De exemplu, dacă N este impar, o permutare și niște schimbări de semne sunt suficiente pentru a transforma relația (4.127.) într-o TF reală de aceeași lungime. Dacă N este par, permutarea următoare:

(4.128.)

transformă relația (4.127.) în:

(4.129.)

Evaluând și simultan, obținem, cu ajutorul relațiilor trigonometrice elementare:

(4.130.)

unde este transformata Fourier a lui (4.128.) și și implică partea reală și respectiv imaginară. Matricea din (4.130.) corespunde unei rotații (sau unei multiplicări complexe).

Evaluarea relației (4.130.) pentru folosește de fiecare dată trei multiplicări și n=N/2 folosește o multiplicare, deci calculul relației (4.130.) are 3N/2-2 multiplicări și 3N/2-3 adunări. Complexitatea unei transformate în cosinus rapide (TCR, de lungime pară este de :

(4.131.)

(4.132.)

Să notăm că acest algoritm nu este optimal în sensul complexității multiplicative. De exemplu, pentru , transformata în cosinus poate fi redusă la un produs polinomial pentru care niște metode optime vor conduce la o diminuare a numărului de multiplicări (dar cu prețul unui număr excesiv de adunări). În general, complexitatea multiplicativă teoretică este legată de cea a transformatei Fourier prin relația următoare:

(4.133.)

și aceasta în virtutea relației (4.129.) care arată că TCR este o parte a TFR de lungime 4N.

Algoritmul pe care l-am descris cu ajutorul relațiilor (4.128.) și (4.130.) este foarte practic și nu folosește decât puține multiplicări în exces ale unui algoritm foarte complex (pentru de exemplu, o multiplicare suplimentară pe când algoritmul optim folosește 29 de adunări suplimentare).

Tabelul 4.5. dă numere de operații pentru transformata în cosinus normalizată (cu factorul pentru ) și asta pentru . Să notăm că o denormalizare prin câștigă două multiplicări (la punctele ) și permite obținerea, de exemplu, TCR de lungime 8 cu ajutorul a numai 11 multiplicări. Aceasta este importantă căci, în practică, TCR de lungimi mici sunt utilizate în codări (tipic, ).

Tabelul 4.5.

Numere de operații pentru transformata în cosinus normalizată de lungime (o versiune denormalizată folosește două multiplicări cel puțin).

În fine, transformata în cosinus bidimensională, definită ca o transformată separabilă cu (4.127.) în fiecare dimensiune, este transformata standard în codarea imaginilor. Dacă numim T matricea transformatei în cosinus, matricea corespunzând transformatei bidimensionale este egală cu:

(4.134.)

similar cu (4.97.). Evident, diverșii algoritmi bidimensionali discutați la secțiunea 4.5 sunt posibili. Derivăm generalizarea algoritmului monodimensional dat prin (4.128.) și (4.130.), căci acest algoritm atinge un număr de operații foarte redus. În notație matriceală, aceste relații indică că T se descompune după cum urmează:

(4.135.)

unde T este o matrice de permutare (4.128.), F este o transformată Fourier reală și D reprezintă rotațiile în (4.131.). Înlocuind (4.135.) în (4.134.) și utilizând proprietatea produsului Kronecker (4.32.):

(4.136.)

Observăm că este calculat cu ajutorul unei permutări bidimensionale , urmată de o TF reală bidimensională , apoi de rotații de ieșire. Aceste rotații transformă patru puncte de ieșire ale TF bidimensionale în patru puncte de ieșire ale TCR bidimensionale.

Dacă notăm și o rotație a unui unghi și respectiv este ușor de verificat că:

(4.137.)

unde

(4.138.)

Aceasta demonstrează că folosește 1,5 multiplicări pe punct, adică rotațiile de ieșire au același cost multiplicativ în una și în două dimensiuni (în fine, rezultatul se generalizează la k dimensiuni). Deci, complexitatea multiplicativă a TCR bidimensională este de:

(4.139.)

Unde cel mai mic provine dintr-un anumit număr de simplificări, ca de exemplu pentru și . Tabelul 4.6 dă numerele de operații pentru calculul TCR bidimensionale, și aceasta pentru algoritmul în linie (coloană ca și pentru algoritmul pe care tocmai l-am derivat).

Tabelul 4.6.

Numărul de operații pentru transformata în cosinus bidimensională de mărime și asta pentru algoritmul în linie/coloană ca și pentru algoritmul direct (versiunea normalizată).

Transformatele în numere întregi

Unul din punctele atrăgătoare ale transformatei Fourier, în afară de utilitatea sa pentru analiza spectrală, este faptul că posedă proprietatea de convoluție. Amintim că la baza unei transformate Fourier de lungime N se găsește un element special care este rădăcina de ordinul N a unității, adică un element de ordin multiplicativ N în corpul numerelor complexe . Ce s-ar întâmpla dacă un element de ordin N există, atunci transformata Fourier poate fi definită, și ea posedă proprietatea de convoluție.

Mai precis, dacă avem un vector de elemente ale GF(q) și că N divide pentru un anume m, atunci elementul de ordin N în este o rădăcină de ordinul N a unității și transformata Fourier de lungime N este definită după cum urmează:

(4.140.)

și transformata inversă este:

(4.141.)

Pentru a găsi inversa lui N, trebuie să îl luăm pe N ca un întreg al corpului GF(q), adică modulo p (unde p este caracteristic GF(q)).

Putem verifica că această transformată Fourier are proprietatea de convoluție, și asta folosind expresiile de mai sus în evaluarea convoluției circulare. Notăm că transformata X(n) a vectorului x(k) trăiește în general într-un corp de extensie, și că deci valorile spectrale sunt legate între ele prin proprietăți de conjugare. Aceasta este deci o generalizare a proprietății de conjugare hermitică a spectrelor Fourier complexe ce provin din secvențe reale.

Sunt două motive principale pentru folosirea transformatelor Fourier pe corpuri finite. Primul este că dacă secvențele de intrare și ieșire pot fi reprezentați exact pe GF(q), atunci convoluția care folosește o transformată Fourier (4.140) și (4.141) definită pe un corp de extensie al GF(q) va fi o soluție exactă, fără orice eroare de rotunjire. Acesta este un avantaj major în raport cu convoluția care folosește transformata Fourier obișnuită (pe ), deoarece atunci rădăcinile de ordin N ale unității sunt aproximate (separat pentru N=2,4) și rezultatul va fi plin de erori de rotunjiri. Dacă secvențele de intrare sunt în virgulă fixă, este ușor să alegem GF(q) suficient de mare pentru a modifica (adapta) atât intrările cât și ieșirile, și este clar că metoda transformatei pe corpuri finite este avantajoasă de vreme ce ea elimină complet zgomotul calculului.

Al doilea motiv pentru anumite lungimi ale cuvintelor și anumite lungimi ale transformatelor, multiplicarea prin rădăcina de ordin N a unității va deveni triviale, de exemplu multiplicarea cu 2. Deci singurele multiplicări restante sunt cele, punct cu punct, cerute pentru convoluție în domeniul transformat, și o convoluție circulară de lungime N se poate calcula cu ajutorul a N multiplicări. Să luăm care nu există decât pentru m luat astfel încât este prim. Aceste numere prime sunt cunoscute sub numele de numerele lui Fermat. O condiție necesară este ca m să fie o putere a lui 2, dar nu este suficientă de vreme ce, de exemplu, este factorizabil. Niște exemple de numere Fermat sunt m=2, 4, 8 și 16 sau 2m +1=5, 17, 257 și respectiv 65537. Există o transformată Fourier pentru toate lungimile care divid 2m , adică n=2, 22 …2m-1 , 2m . Deoarece aceste lungimi sunt puteri ale lui 2, putem utiliza un algoritm Cooley-Tukey, a cărui derivare nu depinde sub nici o formă de corpul ales pentru calcule. Dar, pentru lungimi suficient de mici, o metodă mai interesantă este posibilă. Notăm că elementul 2 are ordinul 2m în și asta căci și deci . De exemplu, în , putem folosi 2 ca o rădăcină de ordinul 32 a unității și putem defini astfel o transformată Fourier de lungime 32.

(4.142.)

O multiplicare printr-o rădăcină a unității este deci o operație de decalaj urmat de o reducere modulo . Această ultimă operație este ea însăți mai degrabă trivială, de vreme ce în și deci adică biții superiori sunt din nou fixați pe 16 poziții și inversați. Transformata dată de relația (4.142) poate fi evaluată fără nici o multiplicare, și o convoluție circulară de lungime 32 nu va utiliza decât 32 multiplicări modulo în domeniul transformat. Bineînțeles, această lungime a transformatei este prea scurtă pentru cea mai mare parte a aplicațiilor, dar o alegere judicioasă a rădăcinii de ordinul N a unității permite adesea reducerea transformatei inițiale la niște mici transformate „triviale”, și asta cu prețul unui număr foarte slab de multiplicări generale.

Un alt caz interesant apare pentru corpul care există dacă este prim. O condiție necesară (dar din nou nu și suficientă) este ca m să fie prim, și de exemplu, este prim pentru m=3, 5, 7, 13, 19 și 31. Aceste numere prime poartă numele de numere Mirsenne. Aritmetica modulo 2m-1 este simplă, de vreme ce, deci raportul este egal cu 1, și biții în exces ai lui 2 sunt doar adunați biților cu pondere slabă. Niște transformate Fourier există pe pentru toate lungimile care divid care nu este deci o putere de 2. De exemplu, este un număr prim, și se factorizează în , ceea ce indică lungimile posibile ale transformatelor Fourier. Notăm că –2 are ordinul 26 pe deoarece . Deci, (-2) permite construirea unei TF de lungime 26.

(4143.)

Evident, multiplicarea prin (-2) este trivială (schimbare a semnului, decalaj, reducția modulo ). Dimpotrivă, 26 nu este o lungime foarte practică.

Notăm că un anumit număr de alte idei, cum ar fi algoritmul Winograd pentru convoluție, și algoritmul Rader pentru TF de o lungime egală cu un număr prim, sunt în mod egal generalizate corpurilor finite. Este suficient să se țină cont de ajustările necesare a faptului că acest corp nu este sau .

Trebuie totuși să remarcăm că aceste metode, deși elegante din puncte de vedere matematic, suferă un număr mare de constrângeri care sunt impuse (pe lungimile cuvintelor și lungimile transformatelor de exemplu). Notăm de asemenea că, înainte de toate, câștigurile sunt interesante în cazul implantărilor materiale (în care aceste aritmetice exotice sunt realizabile cu puțin cost), pe când, pentru implantări în structura logică, avantajele sunt mult mai puțin clare deoarece unitățile aritmetice sunt impuse.

4.11. Algoritmii rapizi pentru soluția sistemelor de ecuații Toeplitz

Matricile Toeplitz

Soluția sistemelor de n ecuații necunoscute este o problemă binecunoscută, și o metodă standard ca eliminarea Gauss folosește operații. În tratarea numerică a semnalului, matricea de dimensiune care caracterizează sistemul de ecuații are adesea o formă specială numită Toeplitz, adică ea este constantă pe lungul diagonalelor. În acest caz, metode mai eficace utilizând a se vedea operații, sunt posibile și constituie subiectul prezentei secțiuni.

Un sistem de ecuații Toeplitz este sub forma:

(4.144.)

sau

(4.145.)

Notăm că matricea Toeplitz A are proprietatea următoare (este persimetrică):

(4.146.)

unde Y este matricea antidiagonală. Adesea, matricea A este simetrică și vectorul Y este de forma:

(4.147.)

adică (4.144.) se scrie astfel:

(4.148.)

În acest caz, soluția (4.148) corespunde rezolvării ecuațiilor Yub-Waken

(4.149.)

care apar, între altele, în problema predicției liniare a parolei utilizate în codaj și în sinteză vocală. În ceea ce urmează, derivăm algoritmii eficace pentru soluția relațiilor (4.144) și (4.148), precum și formulele pentru inversa unei matrici Toeplitz . Notăm de asemenea că soluția sistemelor de ecuații Hankel (constant de-a lungul antidiagonalei) este o problemă centrală în coduri corectoare de erori (pentru codurile Bose-Chaudhum-Hocquenghem). Evident, cum matricile Toeplitz și Hankel sunt legate printr.o permutare (prin Y), toată soluția eficace pentru una din probleme se poate imediat aplica și celeilalte. Metoda clasică în contextul codurilor corectoare de erori este algoritmul Berlekamp-Massey care dă loc unei complexități de ordinul . Acest algoritm poate fi interpretat ca o soluție pentru problema sintezei filtrului autogresiv.

5. Filtrele Gabor

Filtrele Gabor au beneficiat de o atenție considerabilă datorită faptului că, caracteristicile anumitor celule din cortexul vizual al unor mamifere pot fi aproximate (funcționarea acestor celule) prin aceste filtre. Pe lângă aceasta, aceste filtre au arătat că au proprietăți optime de localizare în domeniul spațiu și respectiv frecvență, și ca urmare sunt foarte potrivite pentru problemele de segmentare a structurilor.

Filtrele Gabor sunt utilizate în multe aplicații, ca segmentarea structurilor, detecția unor ținte, managementul dimensiunilor fractale, analiza documentară, detecție, identificare de retină, codări de imagini și reprezentări de imagini.

Un filtru Gabor poate fi vizualizat ca un plan sinusoidal al frecvenței și orientării, modulat de anvelopă gaussiană.

Poate fi scris ca:

unde este o sinusoidă complexă (purtătoarea)

este o funcție bidimensională gaussiană (anvelopa)

Sinusoida complexă e definită după cum urmează

Funcția gaussiană bidimensională este definită de:

Astfel, filtrul Gabor bidimensional (2D) poate fi scris astfel:

Răspunsul în frecvență a filtrului este:

unde .

Acest fapt este echivalent cu translarea funcției gaussiene prin în domeniul frecvență. De aceea funcția Gabor se spune că este o funcție Gauss în domeniul frecvență în poziția și la distanța de de origine și cu o orientare de .

În ecuațiile anterioare, și sunt frecvențele centrale al filtrului Gabor spațial.

Parametrii sunt derivațiile standard ale anvelopei gaussiene pe direcțiile x și y și determinate de banda filtrului.

Fig.5.1. Planul răspunsului în frecvență a filtrului Gabor pentru valori diferite ale , corespunzătoare celor 4 orientări: 0, 45, 90 și 135 grade.

Planul răspunsului în frecvență a filtrului Gabor pentru diferite valori ale lui și , corespund orientărilor 0, 45, 90, 135 (grade) – figura 5.1.

i

Fig.5.2. Ieșirea filtrului Gabor a modelului liniilor simple corespunzătoare celor 4 orientări – 0, 45, 90, 135.

Figura 5.2. arată rezultatul de la ieșirea filtrului Gabor pe 4 direcții, când intrarea este o imagine ce conține linii la aceleași frecvențe doar pentru unghiuri diferite.

Se poate vedea că filtrul cu orientarea are un răspuns corespunzător regiunii de variație când acesta este de asemenea la un unghi . În acest caz, o operație simplă de netezire urmată de un prag sunt suficiente pentru a segmenta imaginea în patru regiuni corespunzătoare liniilor cu cele patru orientări.

Trecând imaginea printr-un filtru Gabor definit prin parametrii (), vom obține toate componentele imaginii ce au energiile concentrate lângă punctul spațial al frecvenței ().

6. Aplicație – Filtrare Gabor recursivă

În această aplicație prezint un algoritm recursiv pentru filtrul Gabor, care realizează prin intermediul unei constante multiplicative o implementare foarte rapidă.

Pentru un semnal compus din N eșantioane, implementarea necesită operatori MADD (Multiply And Add = Multiplică și Adună), al căror număr de operații (operatori) pe fiecare intrare este constant.

Mai mult decât atât, complexitatea este independentă de valorile lui și din nucleul Gabor, și coeficienții ecuației recursive au o soluție simplă și apropiată ca formă, date fiind și .

Implementarea noastră admite nu numai un filtru Gabor de tip „forward” (înainte), dar și un filtru invers care este de asemenea de complexitate (ordin) .

De ce filtrare Gabor?

În timp ce operatorii care pun accent pe informația globală sunt esențiali pentru descrierea unor sisteme fizice dintre cele mai variate, operatorii care privesc informația locală sunt esențiali în analiza semnalelor fizice.

De exemplu sistemele liniare, invariante în timp sunt de obicei analize cu ajutorul transformatelor Fourier sau Laplace, care sunt operații (procese) globale, dar așa cum arată exemplele din figura 1 informația într-un semnal este în general locală.

Fig.1.

Deși sunt multe moduri de abordare a procesărilor de semnale, din punct de vedere al examinării caracterului local al structurii, filtrarea Gabor are câteva proprietăți elegante care o fac foarte potrivită pentru scopul mai sus menționat.

Filtrul Gabor complex, 1-dimensional (unidimensional) – 1D este dat de relația (1).

(1)

Este evident că filtrul Gabor este modulația unui nucleu Gabor prin termenul complex . Utilizarea acestui termen complex semnifică faptul că vom folosi aritmetica complexă în toate calculele.

Pentru a examina structurile locale, vom căuta un filtru a cărui lățime în timp (și spațiu) este mică. Pentru a obține o bună rezoluție a frecvenței, vom căuta un filtru a cărui lățime în frecvență să fie de asemenea mică (îngustă). Anvelopa gaussiană (conturul gaussian) în (1) realizează produsul cu banda de timp cea mai mică posibilă și astfel ne permite să efectuăm o analiză spectrală „locală”.

Din punct de vedere al sincronizării în timp și a selectivității frecvenței, filtrul Gabor este din motive mai sus menționate alegerea optimă.

Mai mult, un număr de autori au evidențiat relația strânsă (similitudinea) între procesarea neurofiziologică a stimulilor tactili și vizuali, și familiile 1D și 2D a filtrelor Gabor. De exemplu, Baker și alții, referindu-se la observațiile privind oscilațiile de tip Gabor în măsurătorile neurofiziologice scrie că „aceste oscilații sunt o latură a sistemului nervos”.

Pe înțelesul nostru, cei mai rapizi algoritmi care au fost descriși sunt în , și pentru un semnal format din N eșantioane, spectrul timp-frecvență are complexitate . Implementarea noastră este a unui filtru Gabor și este bazată pe implementarea recursivă Gaussiană. Din cauza aceasta, vom face un sumar al aspectelor remarcabile ale procesului gaussian.

Analiză a funcției recursive Gauss

Într-o serie de lucrări ale unor cercetători a fost analizată și dezvoltată o metodă de implementare a convoluției de semnale cu nuclee Gauss și utilizarea a două filtre recursive de forma directă

– forward (n, crescător) – ilustrat prin relația (2):

(2) și

– backward – inversă (n descrescător) – ilustrat prin relația (3):

(3)

De notat că în această implementare, ilustrată în figura 2, filtrul „forward” are un răspuns infinit la impuls (IIR), pe care-l putem nota , și filtrul „backward” are un răspuns infinit la impuls .

Concatenarea (legarea) celor două filtre, date de formulele (2) și (3) conduce la un filtru total a cărui transformate Fourier sunt date de:

Cu alte cuvinte, filtrul rezultat are faza „0” (zero)

Fig.2.

Vom aborda filtrarea gaussiană pe un domeniu de frecvență complex și vom utiliza o aproximare rațională a spectrului Gauss:

(4)

unde (5)

De notat că folosim notația de „q” în locul lui „” care este un aspect ce va fi discutat mai târziu.

Polinomul de gradul 6 din denominatorul lui are șase rădăcini, care, datorită structurii speciale a polinomului poate fi scris ca în relația (7) în care

(6)

(7)

Cele trei valori ale lui sunt centrate pe filtrul Gabor. Aceste valori sunt determinate de valorile lui și dacă se alege o altă valoare a lui – pentru a optimiza, de exemplu potrivirea funcției Gauss cu un alt criteriu – atunci valorile lui m se vor modifica de asemenea.

Expresia este scrisă ca un produs a doi termeni;

cu poli în jumătatea stângă a planului, cu poli în jumătatea dreaptă a planului.

Se va utiliza diferența tehnică „backward” (înapoi) și diferența tehnică „forward” (înainte) pentru a transforma cele două filtre analogice stabile caracterizate de și în două filtre digitale stabile și .

Tehnica „backward” (înapoi) este utilizată pentru a transforma în , și tehnica „forward” (înainte) este utilizată pentru a transforma în .

Aceste filtre au transformate Z date de

(8)

(9)

Prelucrarea algebrică folosind „Mathematica” (program) conduce la:

(10)

Cei cinci coeficienți sunt reali și sunt funcții de patru numere date de relațiile din (10).

Primele trei din acestea sunt date de (6).

Ultimul număr q reprezintă valoarea ce trebuie folosită în relațiile (8) și (9) pentru a realiza deviația standard .

Folosind o formulă analitică nouă, care are o îmbunătățire față de formula empirică, relația de legătură între q și pentru filtrul recursiv gaussian este dată de:

(11)

Pentru , anvelopa gaussaină din formula (1) este prea îngustă și subeșantionată și în general se evită.

Algoritmul Gabor

A. De la Gauss la Gabor

Filtrul nostru recursiv utilizează o combinație din implementările noastre Gauss recursive plus o proprietate binecunoscută care face legătura între modulația cu frecvența în domeniul (discret) de timp, și trecerea în domeniul frecvență.

Pentru a fi preciși, dacă semnalul are transformata în Z H(Z), atunci vom avea:

(12)

Înlocuirea lui Z cu reprezintă o rotație a unghiului în jurul punctului Z=(0,0) în planul Z complex.

Consecința directă a acestui fapt este că perechea și din filtrul gauss devin perechea și a filtrului Gabor.

(13)

(14)

unde și sunt polinoame în Z.

Funcția de transfer a filtrului Gabor este dată de formula:

(15)

Coeficientul a fiecărui termen din este chiar conjugata complexă a termenului în al (vezi (13))

Perechea rezultată a ecuației recursive Gabor este:

– „forward”(înainte) – direct – n crescătoare:

(16)

– „backward (înapoi) – invers – n descrescătoare:

(17)

unde

(18)

B. Chestiuni de implementare

Expresia q, descrisă de relația (11), a fost derivată pentru filtrul gaussian.

Pentru filtrul Gabor, problema este mult mai complexă.

O soluție simplă și apropiată de formula (11) nu pare a fi posibilă.

Se va face o analiză regresivă asupra pentru a determina relația dintre ele. Această analiză va conduce la rezultatul empiric:

(19)

Mai departe vom presupune că semnalul unidimensional (1-D) care trebuie filtrat este de lungime N, care este n=1,2,…N.

În direcția cauzală directă („forward”), este important să utilizăm valori proprii pentru primele trei valori din relația (16). Acestea sunt date de:

(20)

Recurența începe la n=4 și merge până la n=N. În direcția inversă (înapoi – „backward”) (anticazuală) este la fel de important să folosim valori proprii pentru „primele” trei valori din (17). Acestea sunt date de

(21)

Recurența care începe la n=N-3 continuă până la n=1. Aplicând acest algoritm unei intrări sub forma unui impuls simplu pentru realizarea echivalentei discrete în timp a relației (1) se obține un răspuns la impuls cum se arată în figura 3 pentru câteva valori ale lui și (figura 3).

Fig.3. Răspunsul complex la impuls g[n] și spectrul Fourier pentru diferite valori ale lui și . Lungimea semnalului N=100. Anvelopele gaussiene au fost trasate pentru comparare.

Ultimul exemplu din figura 3 ne arată un efect interesant. Lățimea totală a semnalului este N=100, și este 20.0. Pentru că răspunsul la impuls să fie suficient de mic la limita semnalului, se impune o scădere a Gaussianului cu aproape 3. Pentru cauze care vor fi clarificate puțin mai târziu, va fi necesar ca răspunsul Gaussian la impuls să fie scăzut cu .

Acest lucru se va produce la ambele capete ale semnalului, care vor fi n=1 și n=N.

Aceasta înseamnă că și N trebuie să îndeplinească condiția . Această condiție nu este satisfăcută în exemplul în care , și acest lucru explică trunchierile capetelor Gaussianului.

Extrapolând acest argument în domeniul frecvență, va fi necesar să se evite distorsionarea spectrului de frecvență, va trebui să avem o valoare astfel încât spectrul gaussian va fi de asemenea scăzută cu aceeași valoare la frecvența maximă . Acest lucru implică faptul că , ceea ce conduce la condiția dată de relația (11). Din aceste două condiții, și anume evitarea aliasing-ului spațial și a celui spectral – conduce la:

(22)

C. Precizia aproximării

Cu toate că figura 3 și relația (22) dau o indicație asupra faptului că modul de alegere a lui poate să producă distorsionări vizibile a rezultatului final, este de asemenea importantă examinarea erorii ce apare chiar în timpul aproximării recursive Gabor.

Pentru a realiza acest lucru (examinarea erorii), vom folosi măsurătorile erorilor care au fost utilizate pentru evaluarea aproximării gaussiene recursive: eroarea absolută și eroarea rădăcină- pătratică:

(23)

unde este diferența complexă dintre răspunsul ideal la impuls Gabor și rezultatul aproximării recursive din formula (17).

Astfel,

Dar

(24)

unde și sunt filtre reale și respectiv Gauss, iar

(25)

unde și sunt aproximările noastre recursive Gabor și respectiv Gauss.

Faptul că acest ultim rezultat dat de relația (25) poate fi factorizat este o consecință directă a structurii noastre de aproximare Gabor dată de relația (12).

Utilizând aceste două rezultate vom vedea că:

(26)

Ultimul termen este eroarea asociată cu aproximarea recursivă gaussiană, o eroare care a fost studiată.

Evaluarea acestei erori măsoară o gamă cuprinsă între așa cum arată în figura 4, rezultatul fiind independent de .

Fig.4. Eroarea măsurată ca o funcție de și . Cele două curbe sunt independente de

D. Filtrare inversă

Procedeul nostru admite de asemenea o implementare directă a filtrului invers. Dacă transformata directă (forward) dată de relația (16) și (17) este dată (în domeniul z) de ecuația:

(27)

atunci filtrul invers este dat de:

(28)

Filtrele și sunt filtre de tip „toți poli” ceea ce înseamnă că inversele acestora sunt filtre de tip „toți-zero” care sunt nerecursive și stabile.

Ecuația va fi dată de:

– forward (direct) – n crescător:

(29)

– backward (invers) – n descrescător:

(30)

E. Complexitatea calculului.

Utilizarea acestei implementări Gabor recursive devine un mod de utilizare a următoarelor „etape”:

Alegerea lui și bazată pe scopul dorit al filtrării;

Utilizarea relației (19) pentru determinarea lui ;

Utilizarea relațiilor (10) și (18) pentru a determina coeficienții;

Utilizarea relațiilor (20) și (21) pentru a inițializa procesul (procedura);

Implementarea filtrului direct (forward) cu relația (16);

Implementarea filtrului inverși (backward) cu relația (17).

Cum (16) și (17) necesită un număr de șapte implementări complexe și șase adunări complexe pentru fiecare punct de ieșire, complexitatea algoritmului de filtrare Gabor pentru N puncte date de intrare este , unde K este o constantă asociată unui MADD complex.

Complexitatea punctelor de ieșire (ieșirilor) este astfel constantă și egală cu .

Aceasta înseamnă că în cadrul acestei constante multiplicative aceasta este cea mai rapidă implementare posibilă, în care fiecare intrare trebuie să fie parcursă cel puțin o dată. Dacă am fi folosit o implementare bazată pe un FFT pentru filtrul Gabor, care presupune faptul că N este un număr compus, atunci complexitatea pe intrare ar fi . În timp ce aceasta poate fi mai mică decât , acest lucru se va întâmpla numai pentru valori mici ale lui N.

De exemplu pentru N=24 =16, filtrul Gabor cu procedeu recursivă este deja mai rapid decât procedeu FFT.

F. Filtrare 2D (bidimensională)

Filtrarea Gabor este utilizată frecvent în procesarea semnalelor multidimensionale. Aplicând procedeul recursiv Gabor descris anterior pentru fiecare dimensiune a semnalului, este posibilă implementarea clasei de filtre Gabor „separabile” (cu separare).

Vom arăta în figura 5 efectul acesta prin trei imagini ce au fost supuse filtrării Gabor aplicate pe câteva direcții.

Filtrul Gabor bidimensional (2D) are un nucleu definit prin:

(31)

Fig.5.

În figura 5, vom avea

Acest filtru poate fi aplicat pe diverse orientări ale lui , precum și pentru valori ale frecvenței (valori) ca și.

Prin figura 5 vom ilustra rezultatele filtrării Gabor pentru valori particular alese ale lui și pentru imaginile ilustrate în figura 1.

În ultimele două rânduri ale figurii 5, am impus condiția . Relația inversă între și e bazată pe proprietăți neurofiziologice ale cortexului vizual.

Alegerea exactă a lui K nu e dificilă și pentru exemplele din figura 5 am ales .

Acest lucru ne duce la o reprezentare considerabilă pentru respectarea relației (22).

Pentru imaginile rezultate a fost aplicată o scalare automatică a contrastului.

Forme de undă Gabor

Este necesară o mențiune pentru ca tehnica descrisă anterior să poată fi folosită pentru reprezentarea formelor de undă ale semnalului Gabor.

O restricție pentru ca o familie de semnale să fie admisă ca forme de undă este aceea că trebuie să aibă media 0, aceasta fiind o componentă dc a valorii 0.

Din rezultatele ilustrate în figura 3, precum și din formula (1) este evident că deși parte imaginară a răspunsului la impuls Gabor are media 0 datorită faptului că este o funcție impară, parte reală a lui este egală și nu are media egală cu 0.

Următorul rezultat care este binecunoscut ne furnizează un mod de a determina valoarea componentei de dc.

(32)

Acest rezultat este destul de complicat dar va fi dat aici pentru ca totul să fie complet:

(33)

unde

(34)

Pentru o valoare a lui m specificată și pentru anumite valori ale lui și , valoarea lui poate fi determinată.

Pentru exemplul dat înainte de figura 3, unde și , .

Extrăgând această valoare din funcția totală de transfer , vom crea o nouă funcție de transfer (o formă de undă Gabor), a cărei valoare este 0 pentru z=1. Această funcție de transfer va avea poli și zerouri:

(35)

Implementarea relației (35) ca un filtru recursiv necesită o abordare alternativă.

Fiecare din polinoamele și sunt polinoame de gradul 3, și de aceea, vor avea trei rădăcini complexe.

Expresia poate fi rescrisă din această cauză ca o operație paralelă în loc de serie, cum se ilustrează în figura 2:

(36)

unde sunt rădăcini ale lui

și sunt rădăcini ale lui .

sunt constante care trebuie determinate prin tehnici standard de operații cu fracții.

Extragerea lui se face prin asocierea cu unul din cei șase termeni, spre exemplu

(37)

Filtrarea formelor de undă din (35) poate fi implementată prin șase operatori paraleli ce vor fi adunați.

Dacă sarcina noastră e analiza formelor de undă, atunci prezența zerourilor nu va fi o problemă, dar dacă intenționăm să recuperăm semnalul original din reprezentarea formelor de undă, atunci tehnica inversă ce implică filtrarea recursivă descrisă anterior nu va funcționa datorită instabilității lui . Vom cerceta un procedeu recursiv pentru recuperarea semnalului din formele de undă (reprezentarea lor) pentru a evita problema aceasta.

7. Concluzii.

Această abordare a filtrării Gabor descrisă aici furnizează un algoritm precis și rapid care este potrivit când valorile lui și sunt cunoscute. Când gama lui trebuie să fie inspectată (de exemplu în spectrul Gabor timp-frecvență), atunci abordarea generală, va trebui reconsiderată.

Complexitatea calcului pentru N intrări este , care este mai rapidă decât o abordare cu un FFt.

În fiecare caz, evitarea erorii alias necesită o condiție de alegere a lui pentru un semnal dat de lungime N, dată de relația (22): .

Algoritmul poate fi aplicat semnalelor multidimensionale când unul singur prezintă interes din clasa filtrelor separabile, ca

(31)

Acest fapt limitează domeniul de aplicabilitate, cu toate acestea există suficiente exemple ce se încadrează în această categorie.

Exemplele includ detectarea liniilor și limitelor obiectelor artificiale.

S-a dezvoltat un filtru Gabor invers, precum și un algoritm pentru reprezentarea formelor de undă.

Încă nu avem un mecanism de producere a recursivităților bazat pe cele două și .

Aceasta înseamnă că semnalul trebuie filtrat pentru fiecare valoare în particular ale pentru a produce rezultatul dorit.

În fiecare domeniu de aplicație, când o relație de genul este cunoscută, o evaluare a tuturor valorilor nu este necesară când se realizează o căutare a lui într-un domeniu dorit.

BIBLIOGRAFIE

Murat, Kunt, Maurice Belanger, Frédéric de Coulon, Claude Gueguen, Techniques modernes du traitement numérique des signaux, Proceedings of IEEE, Paris, France, 1994

Martin Hasler, Gabor Filters, Research Report 1998-2002, Huntsville, Alabama University, May 2002

Ian, T. Young, Lucas J. van Vliet, Michael Ginkel, Recursive Gabor Filtering, IEEE Transactions on Signal Processing, vol.50, No.11, November 2002