Calculer la transformée de Fourier discrète

9

Implémentez la transformée de Fourier discrète (DFT) pour une séquence de n'importe quelle longueur. Cela peut être implémenté sous forme de fonction ou de programme et la séquence peut être donnée soit en argument, soit en utilisant une entrée standard.

L'algorithme calculera un résultat basé sur la DFT standard dans le sens direct. La séquence d'entrée a une longueur Net se compose de [x(0), x(1), ..., x(N-1)]. La séquence de sortie aura la même longueur et se composera de l' [X(0), X(1), ..., X(N-1)]endroit où chacun X(k)est défini par la relation ci-dessous.

DFT

Règles

  • C'est le donc la solution la plus courte l'emporte.
  • Les fonctions intégrées qui calculent la DFT dans les directions avant ou arrière (également appelées inverses) ne sont pas autorisées.
  • Les inexactitudes en virgule flottante ne seront pas prises en compte contre vous.

Cas de test

DFT([1, 1, 1, 1]) = [4, 0, 0, 0]
DFT([1, 0, 2, 0, 3, 0, 4, 0]) = [10, -2+2j, -2, -2-2j, 10, -2+2j, -2, -2-2j]
DFT([1, 2, 3, 4, 5]) = [15, -2.5+3.44j, -2.5+0.81j, -2.5-0.81j, -2.5-3.44j]
DFT([5-3.28571j, -0.816474-0.837162j, 0.523306-0.303902j, 0.806172-3.69346j, -4.41953+2.59494j, -0.360252+2.59411j, 1.26678+2.93119j] = [2, -3j, 5, -7j, 11, -13j, 17]

Aidez-moi

Il y avait un défi précédent pour trouver la DFT en utilisant un algorithme FFT pour des séquences de longueurs égales à une puissance de 2. Vous pouvez y trouver quelques astuces qui pourraient vous aider ici. Gardez à l'esprit que ce défi ne vous limite à aucune complexité et nécessite également que votre solution fonctionne pour des séquences de n'importe quelle longueur.

miles
la source

Réponses:

2

Gelée , 16 15 octets

LR’µ×'÷L-*²³÷S€

Essayez-le en ligne!

Comment ça fonctionne

LR’µ×'÷L-*²³÷S€  Main link. Argument [x(0), ..., x(N-1)].

L                Length; yield N.
 R               Range; yield [1, ..., N].
  ’              Decrement; yield [0, ..., N-1].
   µ             Begin a new, monadic chain. Argument: [0, ..., N-1]
    ×'           Spawned multiply [0, ..., N-1] with itself, yielding the matrix
                 of all possible products k×n.
      ÷L         Divide each product by N.
        -*       Compute (-1)**(kn÷L) for each kn.
          ²      Square each result, computing (-1)**(2kn÷L).
           ³÷    Divide [x(0), ..., x(N-1)] by the results.
             S€  Compute the sum for each row, i.e., each X(k).
Dennis
la source
ninja'd :(
Leaky Nun
5

Python 3, 77 octets

lambda x,e=enumerate:[sum(t/1j**(4*k*n/len(x))for n,t in e(x))for k,_ in e(x)]

Testez-le sur Ideone .

Le code utilise la formule équivalente

formule

Dennis
la source
Wow, des chiffres énormes. C'est agréable de voir les formules équivalentes qui peuvent permettre un code plus court.
miles
4

J, 30 20 octets

3 octets grâce à @miles .

Utilise le fait que e^ipi = -1.

La formule devient X_k = sum(x_n / ((-1)^(2nk/N))).

+/@:%_1^2**/~@i.@#%#

Usage

>> DCT =: +/@:%_1^2**/~@i.@#%#
>> DCT 1 2 3 4 5
<< 15 _2.5j3.44095 _2.5j0.812299 _2.5j_0.812299 _2.5j_3.44095

>>est STDIN et <<STDOUT.

"Les inexactitudes en virgule flottante ne seront pas prises en compte contre vous."

Leaky Nun
la source
3

MATL , 20 16 octets

-1yn:qt!Gn/E*^/s

L'entrée est un vecteur colonne, c'est-à-dire remplacer les virgules par des points-virgules:

[1; 1; 1; 1]
[1; 0; 2; 0; 3; 0; 4; 0]
[1; 2; 3; 4; 5]
[5-3.28571j; -0.816474-0.837162j; 0.523306-0.303902j; 0.806172-3.69346j; -4.41953+2.59494j; -0.360252+2.59411j; 1.26678+2.93119j] 

Cela utilise la formule dans la réponse de Leaky Nun , basée sur les faits que exp ( ) = −1, et que l'opératon de puissance de MATL avec un exposant non entier produit (comme dans la plupart des langages de programmation) le résultat de la branche principale .

Essayez-le en ligne!

En raison de l'espacement étrange d'Octave avec des nombres complexes, les parties réelles et imaginaires sont séparées par un espace, tout comme les différentes entrées du vecteur résultant. Si cela semble trop moche, ajoutez un !à la fin ( 17 octets ) pour avoir chaque entrée de la sortie dans une ligne différente.

Explication

-1      % Push -1
y       % Get input implicitly. Push a copy below and one on top of -1
n:q     % Row vector [0 1 ... N-1] where N is implicit input length
t!      % Duplicate and transpose: column vector
Gn      % Push input length
/       % Divide
E       % Multiply by 2
*       % Multiply, element-wise with broadcast. Gives the exponent as a matrix
^       % Power (base is -1), element-wise. Gives a matrix
/       % Divide matrix by input (column vector), element-wise with broadcast
s       % Sum
Luis Mendo
la source
2

Pyth, 30

ms.e*b^.n1****c_2lQk.n0d.j0)QU

Suite de tests

Approche très naïve, juste une implémentation de la formule. Se heurte à divers problèmes mineurs de virgule flottante avec des valeurs qui devraient être des inverses additifs s'ajoutant pour aboutir à des valeurs légèrement différentes de zéro.

Curieusement, .jcela ne semble pas fonctionner sans argument, mais je ne sais pas si je l'utilise correctement.

FryAmTheEggman
la source
1
Félicitations pour 10k !!
Luis Mendo
2

Pyth, 18 octets

Utilise le fait que e^ipi = -1.

La formule devient X_k = sum(x_n / ((-1)^(2nk/N))).

ms.ecb^_1*cyklQdQU

Suite de tests.

Leaky Nun
la source
2

Python 2, 78 octets

l=input();p=1
for _ in l:print reduce(lambda a,b:a*p+b,l)*p;p*=1j**(4./len(l))

Le polynôme est évalué pour chaque puissance pde 1j**(4./len(l)).

L'expression reduce(lambda a,b:a*p+b,l)évalue le polynôme donné par lsur la valeur pvia la méthode de Horner:

reduce(lambda a,b:a*10+b,[1,2,3,4,5])
=> 12345

Sauf que la liste d'entrée est inversée, avec un terme constant à la fin. Nous pourrions l'inverser, mais parce que p**len(l)==1pour les coefficients de Fourier, nous pouvons utiliser un hack d'inversionp et multiplier le résultat entier par p.

Quelques tentatives de longueur égale:

l=input();i=0
for _ in l:print reduce(lambda a,b:(a+b)*1j**i,l,0);i+=4./len(l)

l=input();i=0
for _ in l:print reduce(lambda a,b:a*1j**i+b,l+[0]);i+=4./len(l)

En fonction pour 1 octet de plus (79):

lambda l:[reduce(lambda a,b:a*1j**(i*4./len(l))+b,l+[0])for i in range(len(l))]

Une tentative de récursivité (80):

f=lambda l,i=0:l[i:]and[reduce(lambda a,b:(a+b)*1j**(i*4./len(l)),l,0)]+f(l,i+1)

Simulation itérative reduce(80):

l=input();p=1;N=len(l)
exec"s=0\nfor x in l:s=s*p+x\nprint s*p;p*=1j**(4./N);"*N
xnor
la source
2

C (gcc) , 86 78 octets

k;d(a,b,c)_Complex*a,*b;{for(k=c*c;k--;)b[k/c]+=a[k%c]/cpow(1i,k/c*k%c*4./c);}

Essayez-le en ligne!

Cela suppose que le vecteur de sortie est remis à zéro avant utilisation.

plafond
la source
1

Python 2, 89 octets

Utilise le fait que e^ipi = -1.

La formule devient X_k = sum(x_n / ((-1)^(2nk/N))).

lambda a:[sum(a[x]/(-1+0j)**(x*y*2./len(a))for x in range(len(a)))for y in range(len(a))]

Ideone it!

Leaky Nun
la source
Publiez cela comme une réponse séparée!
Leaky Nun
D'accord, si vous le dites.
Dennis
1

Mathematica, 49 48 47 octets

Total[I^Array[4(+##-##-1)/n&,{n=Length@#,n}]#]&

Basé sur la formule de la solution @Dennis .

miles
la source
1

Axiome, 81 octets

g(x)==(#x<2=>x;[reduce(+,[x.j/%i^(4*k*(j-1)/#x)for j in 1..#x])for k in 0..#x-1])

en utilisant la formule que quelqu'un poste ici. Résultats

(6) -> g([1,1,1,1])
   (6)  [4,0,0,0]
                                    Type: List Expression Complex Integer
(7) -> g([1,2,3,4])
   (7)  [10,- 2 + 2%i,- 2,- 2 - 2%i]
                                    Type: List Expression Complex Integer
(8) -> g([1,0,2,0,3,0,4,0])
   (8)  [10,- 2 + 2%i,- 2,- 2 - 2%i,10,- 2 + 2%i,- 2,- 2 - 2%i]
                                    Type: List Expression Complex Integer
(11) -> g([1,2,3,4,5])
   (11)
        5+--+4       5+--+3    5+--+2      5+--+
        \|%i   + 5%i \|%i   - 4\|%i   - 3%i\|%i  + 2
   [15, --------------------------------------------,
                           5+--+4
                           \|%i
    5+--+4       5+--+3    5+--+2      5+--+
    \|%i   + 3%i \|%i   - 5\|%i   - 2%i\|%i  + 4
    --------------------------------------------,
                       5+--+4
                       \|%i
    5+--+4       5+--+3    5+--+2      5+--+
    \|%i   + 4%i \|%i   - 2\|%i   - 5%i\|%i  + 3
    --------------------------------------------,
                       5+--+4
                       \|%i
    5+--+4       5+--+3    5+--+2      5+--+
    \|%i   + 2%i \|%i   - 3\|%i   - 4%i\|%i  + 5
    --------------------------------------------]
                       5+--+4
                       \|%i
                                    Type: List Expression Complex Integer
(12) -> g([1,2,3,4,5.])
   (12)
   [15.0, - 2.5 + 3.4409548011 779338455 %i, - 2.5 + 0.8122992405 822658154 %i,
    - 2.5 - 0.8122992405 822658154 %i, - 2.5 - 3.4409548011 779338455 %i]
                                      Type: List Expression Complex Float
RosLuP
la source