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

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

ωN=e

2iπ
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

   pjexp(ijx),p(xk)=fk

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

   N−1
   ∑
   j=0

   pjexp(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

   cnexp(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)