#include<stdio.h>
#include<stdlib.h>
#include<malloc.h>
#include<math.h>
#include"complexe.h"
#include"fourier.h"
#include"linkedList.h"
/*------------------------------------------*/
/*
Retourne les racines n-eme de c, O(n)
*/
linkedList* racinesComplexe(complexe c, unsigned long n)
{
linkedList* l = linkedListCreate();
pComplexe p;
unsigned long i;
for(i = 0 ; i < n ; i++)
{
p = (pComplexe)malloc(sizeof(complexe));
if (p == NULL)
exit(1);
*p = makeExpComplexe(pow(c.mod, 1./n), (c.arg/n) + (i)*2*M_PI/n);
linkedListAppend(l, p);
}
return l;
}
/*------------------------------------------*/
/*
Retourne les racines n-eme de 1, O(n)
*/
linkedList* uniteRacinesComplexe(unsigned long n)
{
return racinesComplexe(makeAffComplexe(1.0, 0.0), n);
}
/*------------------------------------------*/
/*
Donne l'image de x par polynome. Les coefficients du polynome sont
ranges par ordre de degre croissant. O(n).
*/
complexe horner(link* polynome, complexe x)
{
if (polynome == NULL)
return makeAffComplexe(0.0, 0.0);
return addComplexe(*(pComplexe)polynome->data,
multComplexe(x,
horner(polynome->next,
x)));
}
/*------------------------------------------*/
/*
Detruit complexe pointe par p.
*/
void deallocateComplexe(void* p)
{
free((pComplexe)p);
}
/*------------------------------------------*/
/*
Retourne la transformee discrete de Fourier du vecteur l,
version iterative en O(n^2)
*/
linkedList* transformeeFourierIt(linkedList* polynome)
{
linkedList* FT = linkedListCreate();
linkedList* W_n = uniteRacinesComplexe(linkedListGetSize(polynome));
void* evaluate(void* w)
{
pComplexe p = (pComplexe)malloc(sizeof(complexe));
if (p == NULL)
exit(0);
*p = horner(linkedListGetFirst(polynome),
*(pComplexe)w);
return p;
}
FT = linkedListMap(W_n, evaluate);
linkedListDestroy(W_n, deallocateComplexe);
return FT;
}
/*------------------------------------------*/
/*
Ne fait rien...
*/
void doNothing(void* w)
{
}
/*------------------------------------------*/
/*
Retourne la transformee discrete de Fourier du vecteur l,
version recursive en O(n.log n)
*/
linkedList* FFT(linkedList* polynome)
{
// creation des listes
linkedList* FT = linkedListCreate(); // resultat
linkedList* odd = linkedListCreate(); // coefficients d'indice pair
linkedList* even = linkedListCreate();// coefficients d'indice impair
linkedList* a_0 = linkedListCreate(); // TF de odd
linkedList* a_1 = linkedListCreate(); // TF de even
link* l; // pour parcourir FT
complexe w = makeAffComplexe(1.0, 0.0);
complexe c = makeExpComplexe(1.0, 2*M_PI/linkedListGetSize(polynome));
int isEven = 1;
// cas de base
if (linkedListGetSize(polynome) == 1)
{
linkedListAppend(FT,
copyComplexe(*(pComplexe)
(linkedListGetFirst(polynome)->data)));
return FT;
}
// autres cas
void split(void* coeff)
{
if (isEven)
linkedListAppend(even, coeff);
else
linkedListAppend(odd, coeff);
isEven = !isEven;
}
// separation de polynome
linkedListApply(polynome, split);
// appels recursifs
a_0 = FFT(even);
a_1 = FFT(odd);
void evenValues(void* v)
{
linkedListAppend(FT, copyComplexe(*(pComplexe)v));
}
// ajout dans FT des valeurs de a^[0]
linkedListApply(a_0, evenValues);
linkedListApply(a_0, evenValues);
l = linkedListGetFirst(FT);
void oddValues(void* v)
{
pComplexe ft = (pComplexe)(l->data);
*ft = addComplexe(*ft, multComplexe(w, *(pComplexe)v));
w = multComplexe(c, w);
l = l->next;
}
// ajout dans FT des w_n^p * a^[1]_p
linkedListApply(a_1, oddValues);
linkedListApply(a_1, oddValues);
// liberation de la memoire
linkedListDestroy(odd, doNothing);
linkedListDestroy(even, doNothing);
linkedListDestroy(a_0, deallocateComplexe);
linkedListDestroy(a_1, deallocateComplexe);
return FT;
}
/*------------------------------------------*/
/*
Retourne la transformee discrete de Fourier du vecteur l,
version recursive en O(n.log n)
*/
linkedList* inverseFFT(linkedList* polynome)
{
linkedList* FT = linkedListCreate();
linkedList* odd = linkedListCreate();
linkedList* even = linkedListCreate();
linkedList* a_0 = linkedListCreate();
linkedList* a_1 = linkedListCreate();
link* l;
complexe w = makeAffComplexe(1.0, 0.0);
complexe c = makeExpComplexe(1.0, -2.*M_PI/linkedListGetSize(polynome));
int isEven = 1;
if (linkedListGetSize(polynome) == 1)
{
linkedListAppend(FT,
copyComplexe(*(pComplexe)
(linkedListGetFirst(polynome)->data)));
return FT;
}
void split(void* coeff)
{
if (isEven)
linkedListAppend(even, coeff);
else
linkedListAppend(odd, coeff);
isEven = !isEven;
}
linkedListApply(polynome, split);
a_0 = inverseFFT(even);
a_1 = inverseFFT(odd);
void evenValues(void* v)
{
linkedListAppend(FT, copyComplexe(*(pComplexe)v));
}
linkedListApply(a_0, evenValues);
linkedListApply(a_0, evenValues);
l = linkedListGetFirst(FT);
void oddValues(void* v)
{
pComplexe ft = (pComplexe)(l->data);
complexe value = multComplexe(w, *(pComplexe)v);
value = multReelComplexe(1./2., addComplexe(*ft, value));
*ft = value;
w = multComplexe(c, w);
l = l->next;
}
linkedListApply(a_1, oddValues);
linkedListApply(a_1, oddValues);
void doNothing(void* w){}
linkedListDestroy(odd, doNothing);
linkedListDestroy(even, doNothing);
void deallocateComplexe(void* p)
{
free((pComplexe)p);
}
linkedListDestroy(a_0, deallocateComplexe);
linkedListDestroy(a_1, deallocateComplexe);
return FT;
}
/*------------------------------------------*/
/*
Retourne la transformee discrete de Fourier inverse
du vecteur l, version iterative en O(n^2)
*/
linkedList* transformeeInverseFourierIt(linkedList* polynome)
{
linkedList* FT = linkedListCreate();
int n = linkedListGetSize(polynome);
linkedList* W_n = uniteRacinesComplexe(n);
void conjugue(void* w)
{
*(pComplexe)w = conjugueComplexe(*(pComplexe)w);
}
linkedListApply(W_n, conjugue);
void* evaluate(void* w)
{
pComplexe p = (pComplexe)malloc(sizeof(complexe));
if (p == NULL)
exit(0);
*p = horner(linkedListGetFirst(polynome),
*(pComplexe)w);
return p;
}
FT = linkedListMap(W_n, evaluate);
void divise(void* w)
{
*(pComplexe)w = multReelComplexe(1./n, *(pComplexe)w);
}
linkedListApply(FT, divise);
void deallocateComplexe(void* p)
{
free((pComplexe)p);
}
linkedListDestroy(W_n, deallocateComplexe);
return FT;
}
/*------------------------------------------*/
linkedList* simpleList(int n)
{
if (n == 0)
return linkedListCreate();
linkedList* l = simpleList(n - 1);
linkedListAppend(l, copyComplexe(makeAffComplexe(n, 0.0)));
return l;
}
/*------------------------------------------*/
/*
int main()
{
int k = 512;
void printC(void* v)
{
printAffComplexe(*((pComplexe)v));
printf("\n");
}
printf("liste\n");
linkedListApply(simpleList(k),
printC);
printf("FFT de liste\n");
linkedListApply(FFT(simpleList(k)),
printC);
printf("FFT-1 o FFT de liste\n");
linkedListApply(inverseFFT(FFT(simpleList(k))),
printC);
return 0;
}
*/