### 1.22.2 Διακριτός Μετασχηματισμός `Fourier```

Έστω N ένας ακέραιος. Ο Διακριτός Μετασχηματισμός `Fourier (DFT)` είναι ένας μετασχηματισμός FN που ορίζεται στο σύνολο των περιοδικών ακολουθιών περιόδου N. Εξαρτάται απ’ την επιλογή μιας αρχικής N-οστής ρίζας της μονάδας ωN. Εάν ο `DFT` ορίζεται σε ακολουθίες με μιγαδικούς συντελεστές, έχουμε:

ωN=e

2 i π
N

Εάν x είναι μια περιοδική ακολουθία περιόδου N, που ορίζεται από το διάνυσμα x=[x0,x1,...xN−1] τότε ο FN(x)=y είναι μια περιοδική ακολουθία περιόδου N, που ορίζεται από:

(FN,ωN(x))k=yk=

N−1
∑
j=0

xjωN−k· j, k=0..N−1

όπου ωN είναι μια αρχική N-οστή ρίζα της μονάδας. Mε τον Γρήγορο Μετασχηματισμό `Fourier (FFT)` o διακριτός Μετασχηματισμός `Fourier` μπορεί να υπολογιστεί γρηγορότερα από τον υπολογισμό του καθενός yk ατομικά. Το `Xcas` υλοποιεί τον αλγόριθμο `FFT` για να υπολογίσει τον διακριτό μετασχηματισμό `Fourier` μόνο εάν το N είναι δύναμη του 2.

#### Ιδιότητες του διακριτού μετασχηματισμού `Fourier`

Ο Διακριτός Μετασχηματισμός `Fourier` FN είναι ένας αμφιμονοσήμαντος και επί μετασχηματισμός σε περιοδικές ακολουθίες τέτοιες ώστε :

FN,ωN−1

=

1
N

FN,ωN−1

=

1
N

FN

στο ℂ

δηλαδή :

(FN−1(x))k=

1
N

N−1
∑
j=0

xjωNk· j

Στο `Xcas` ο Διακριτός Μετασχηματισμός `Fourier` και ο αντίστροφός του δηλώνονται με `fft` και `ifft`:

`fft(x)`= FN(x), `ifft(x)`= FN−1(x)

Ορισμοί
Έστω x και y δύο περιοδικές ακολουθίες με περίοδο N.

- Το γινόμενο `Hadamard` (συμβολισμός ·) ορίζεται με:

  (x · y)k = xk yk

- το γινόμενο συνέλιξης (συμβολισμός *) ορίζεται με:

  (x * y)k=

  N−1
  ∑
  j=0

  xjyk−j

Ιδιότητες :

N*FN(x · y) | = | FN(x) * FN(y)
FN(x * y) | = | FN(x) · FN(y)

#### Εφαρμογές

1. Τιμή ενός πολυωνύμου
   Ορίζουμε ένα πολυώνυμο P(x)=∑j=0N−1cjxj με το διάνυσμα των συντελεστών του c:=[c0,c1,..cN−1], όπου μηδενικά μπορούν να προστεθούν έτσι ώστε το N να είναι δύναμη του 2.
   - Υπολογισμός των τιμών του P(x) στα σημεία

     x=ak=ωN−k=exp(

     −2ikπ
     N

     ), k=0..N−1

     Αυτός είναι απλά ο διακριτός μετασηματισμός `Fourier` του c αφού,

     P(ak)=

     N−1
     ∑
     j=0

     cj(ωN−k)j=FN(c)k

     Εισάγετε, για παράδειγμα :
     `P(x):=x+x```^```2; w:=i`
     Εδώ οι συντελεστές του P είναι [0,1,1,0], N=4 και ω=exp(2iπ/4)=i.
     Είσοδος :
     `fft([0,1,1,0])`
     Έξοδος :
     `[2,-1-i,0,-1+i]`
     επομένως,
     - `P(0)=2`,
     - `P(-i)=P(w```^```-1)=-1-i`,
     - `P(-1)=P(w```^```-2)=0`,
     - `P(i)=P(w```^```-3)=-1+i`.
   - Υπολογισμός των τιμών του P(x) στα

     x=bk=ωNk=exp(

     2ikπ
     N

     ), k=0..N−1

     Αυτό είναι N επί τον αντίστροφο μετασχηματισμό `Fourier` του c αφού,

     P(ak)=

     N−1
     ∑
     j=0

     cj(ωNk)j=NFN−1(c)k

     Εισάγετε, για παράδειγμα :
     `P(x):=x+x```^```2` και `w:=i`
     Εδώ, οι συντελεστές του P είναι [0,1,1,0], N=4 και ω=exp(2iπ/4)=i.
     Είσοδος :
     `4*ifft([0,1,1,0])`
     Έξοδος :
     `[2,-1+i,0,-1-i]`
     επομένως :
     - `P(0)=2`,
     - `P(i)=P(w```^```1)=-1+i`,
     - `P(-1)=P(w```^```2)=0`,
     - `P(-i)=P(w```^```3)=-1-i`.
     Ασφαλώς, βρίσκουμε τις ίδιες τιμές όπως και παραπάνω...
2. Τριγωνομετρική παρεμβολή
   Έστω f μια περιοδική συνάρτηση περιόδου 2π, και έστω ότι f(2kπ/N)=fk for k=0..(N−1). Βρείτε ένα τριγωνομετρικό πολυώνυμο p που παρεμβάλει την f στα xk=2kπ/N, δηλαδή βρείτε pj, j=0..N−1 τέτοια ώστε

   p(x)=

   N−1
   ∑
   j=0

   pj exp(ijx), p(xk)=fk

   Αντικαθιστώντας στο p(x) το xk με την τιμή του έχουμε:

   N−1
   ∑
   j=0

   pj exp(i

   j2kπ
   N

   ) = fk

   Μ’ άλλα λόγια, (fk) είναι ο αντίστροφος μετασχηματισμός `Fourier DFT` του (pk), και επομένως

   (pk)=

   1
   N

   FN( (fk) )

   Εάν η συνάστηση f είναι πραγματική, p−k=pk, και επομένως ανάλογα αν το N είναι άρτιο ή περιττό:

   p(x)

   =

   p0+ 2 ℜ(

   N
   2

   −1

   ∑
   k=0

   pkexp(ikx))+ℜ(p

   N
   2

   exp(i

   Nx
   2

   ))

   p(x)

   =

   p0+ 2 ℜ(

   N−1
   2

   ∑
   k=0

   pkexp(ikx))

3. Σειρές `Fourier`
   Έστω f μια περιοδική συνάρτηση περιόδου 2π, τέτοια ώστε :

   f(xk)=yk, xk=

   2kπ
   N

   , k=0..N−1

   Έστω επίσης ότι η σειρά `Fourier` της f συγκλίνει στην f (αυτό θα ισχύει εάν για παράδειγμα f είναι συνεχής). Εάν το N είναι μεγάλο, μια καλή προσέγγιση της f θα δίνεται από:

   ∑

   −

   N
   2

   ≤ n<

   N
   2

   cn exp(inx)

   Επομένως θέλουμε μια αριθμητική προσέγγιση των

   cn=

   1
   2π

   ∫

   2π
   0

   f(t)exp(−int)dt

   Η αριθμητική τιμή του ολοκληρώματος ∫02πf(t)exp(−int)dt μπορεί να υπολογιστεί από τον τραπεζοειδή κανόνα (σημειώνουμε ότι ο αλγόριθμος `Romberg` δεν θα δούλευε εδώ, επειδή το ανάπτυγμα `Euler Mac Laurin` έχει μηδενικούς συντελεστές, αφού η ολοκληρωτέα συνάρτηση είναι περιοδική, και επομένως όλες της οι παράγωγοι έχουν την ίδια τιμή στο 0 και στο 2π). Εάν c_n είναι η αριθμητική τιμή της cn που παίρνουμε από τον τραπεζοειδή κανόνα, τότε

   c_n=

   1
   2π

   2π
   N

   N−1
   ∑
   k=0

   ykexp(−2i

   nkπ
   N

   ), −

   N
   2

   ≤ n<

   N
   2

   Πράγματι, αφού xk=2kπ/N και f(xk)=yk:

   f(xk)exp(−inxk)

   =

   ykexp(−2i

   nkπ
   N

   ),

   f(0)exp(0)=f(2π)exp(−2i

   nNπ
   N

   )

   =

   y0=yN

   Επομένως :

   [c0,..c

   N
   2

   −1

   ,c

   −N
   2

   ,..c−1]=

   1
   N

   FN([y0,y1...y(N−1)])

   <αφού
   - εάν n≥0, cn=yn
   - εάν n<0 cn=yn+N
   - ωN=exp(2iπ/N), τότε ωNn=ωNn+N

   Ιδιότητες

   - Οι συντελεστές του τριγωνομετρικού πολυωνύμου που παρεμβάλει την f στα x=2kπ/N είναι

     pn=cn, −

     N
     2

     ≤ n<

     N
     2

   - Εάν η f είναι ένα τριγωνομετρικό πολυώνυμο P βαθμού m≤ N/2, τότε

     f(t)=P(t)=

     m−1
     ∑
     k=−m

     ckexp(2ikπ t)

     το τριγωνομετρικό πολυώνυμο που παρεμβάλει την f=P είναι το P, η αριθμητική προσέγγιση των συντελεστών είναι στην πραγματικότητα ακριβής (cn=cn).
   - Πιο γενικά, μπορούμε να υπολογίσουμε το cn−cn.
     Υποθέτοντας ότι η f ισούται με την σειρά Fourier της, δηλαδή ότι :

     f(t)=

     +∞
     ∑
     m=−∞

     cmexp(2iπ mt),

     +∞
     ∑
     m=−∞

     |cm|<∞

     Τότε :

     f(xk)=f(

     2kπ
     N

     )=yk=

     +∞
     ∑
     m=−∞

     cmωNkm, c_n=

     1
     N

     N−1
     ∑
     k=0

     ykωN−kn

     Αντικαταστείστε το yk με την τιμή του στο c_n:

     c_n=

     1
     N

     N−1
     ∑
     k=0

     +∞
     ∑
     m=−∞

     cmωNkmωN−kn

     Εάν m≠ n (mod N ), ωNm−n είναι μια N-στη ρίζα της μονάδος διάφορη του 1, και επομένως:

     ωN(m−n)N=1,

     N−1
     ∑
     k=0

     ωN(m−n)k=0

     Άρα, εάν m−n είναι ένα πολλαπλάσιο του N (m=n+l· N) τότε ∑k=0N−1ωNk(m−n)=N, διαφορετικά ∑k=0N−1ωNk(m−n)=0. Αντιστρέφοντας τα δυο αθροίσματα, έχουμε

     c_n

     =

     1
     N

     +∞
     ∑
     m=−∞

     cm

     N−1
     ∑
     k=0

     ωNk(m−n)

     =

     +∞
     ∑
     l=−∞

     c(n+l· N)

     | = | ...cn−2· N+cn−N+cn+cn+N+cn+2· N+.....

     Συμπέρασμα: εάν |n|<N/2, η διαφορά c_n−cn είναι ένα άθροισμα του cj μεγάλων δεικτών (τουλάχιστον N/2 σε απόλυτη τιμή), και θα είναι μικρή (εξαρτάται από την τάξη σύγκλισης των σειρών Fourier).

   Παράδειγμα, εισάγετε

   `f(t):=cos(t)+cos(2*t)`
   `x:=f(2*k*pi/8)$(k=0..7)`

   Έξοδος :

   `x={2,(√``2``)/2,-1,-((√``2``)/2),0,-((√``2``)/2),-1,(√``2``)/2}`
   `fft(x)=[0.0,4.0,4.0,0.0,0.0,0.0,4.0,4.0]`

   Μετά από διαίρεση με N=8, έχουμε :

   c0=0,c1=4.0/8,c2=4.0/2,c3=0.0,
   c−4=0,c−3=0,c−2=4.0/8,=c−1=4.0/8

   Επομένως, bk=0 και ak=c−k+ck ισούται με 1 εάν k=1, 2 και 0 σε κάθε άλλη περίπτωση.

4. Γινόμενο Συνέλιξης
   Εάν P(x)=∑j=0n−1ajxj και Q(x)=∑j=0m−1bjxj δίνονται από τα διανύσματα των συντελεστών τους a=[a0,a1,..an−1] and b=[b0,b1,..bm−1], υπολογίζουμε το γινόμενο αυτών των δύο πολυωνύμων χρησιμοποιώντας τον `DFT`. Το γινόμενο των δύο πολυωνύμων είναι το γινόμενο συνέλιξης της περιοδικής σειράς των συντελεστών τους εάν η περίοδος είναι μεγαλύτερη ή ίση από (n+m). Επομένως συμπληρώνουμε το a (αντίστοιχα το b) με m+p (αντίστοιχα n+p) μηδενικά, όπου p επιλέγεται τέτοιο ώστε το N=n+m+p να είναι δύναμη του 2. Εάν a=[a0,a1,..an−1,0..0] και b=[b0,b1,..bm−1,0..0], τότε:

   P(x)Q(x)=

   n+m−1
   ∑
   j=0

   (a*b)jxj

   Υπολογίζουμε τα FN(a), FN(b), και μετά ab=FN−1(FN(a)· FN(b)) χρησιμοποιώντας τις ιδιότητες

   NFN(x · y)=FN(x) * FN(y), FN(x * y)=FN(x) · FN(y)