Calcule la transformada discreta de Fourier

9

Implemente la Transformada discreta de Fourier (DFT) para una secuencia de cualquier longitud. Esto puede implementarse como una función o un programa y la secuencia se puede dar como un argumento o utilizando una entrada estándar.

El algoritmo calculará un resultado basado en DFT estándar en la dirección hacia adelante. La secuencia de entrada tiene longitud Ny consta de [x(0), x(1), ..., x(N-1)]. La secuencia de salida tendrá la misma longitud y consiste en [X(0), X(1), ..., X(N-1)]donde cada uno X(k)está definido por la siguiente relación.

DFT

Reglas

  • Este es el por lo que gana la solución más corta.
  • No se permiten las unidades integradas que calculan el DFT en direcciones hacia adelante o hacia atrás (también conocido como inverso).
  • Las imprecisiones de punto flotante no se contarán en su contra.

Casos de prueba

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]

Ayuda

Hubo un desafío anterior para encontrar el DFT usando un algoritmo FFT para secuencias con longitudes iguales a una potencia de 2. Puede encontrar algunos trucos allí que podrían ayudarlo aquí. Tenga en cuenta que este desafío no lo limita a ninguna complejidad y también requiere que su solución funcione para secuencias de cualquier longitud.

millas
fuente

Respuestas:

2

Jalea , 16 15 bytes

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

Pruébalo en línea!

Cómo funciona

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
fuente
ninja'd :(
Leaky Nun
5

Python 3, 77 bytes

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

Pruébalo en Ideone .

El código usa la fórmula equivalente

fórmula

Dennis
fuente
Wow, figuras gigantescas. Es agradable ver las fórmulas equivalentes que pueden permitir un código más corto.
millas
4

J, 30 20 bytes

3 bytes gracias a @miles .

Utiliza el hecho de que e^ipi = -1.

La fórmula se convierte X_k = sum(x_n / ((-1)^(2nk/N))).

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

Uso

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

donde >>es STDIN y <<es STDOUT.

"Las imprecisiones de coma flotante no se contarán en su contra".

Monja permeable
fuente
3

MATL , 20 16 bytes

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

La entrada es un vector de columna, es decir, reemplazar comas por punto y coma:

[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] 

Esto usa la fórmula en la respuesta de Leaky Nun , basada en los hechos que exp ( ) = −1, y que la operación de poder de MATL con exponente no entero produce (como en la mayoría de los lenguajes de programación) el resultado de la rama principal .

Pruébalo en línea!

Debido al extraño espacio de Octave con números complejos, las partes real e imaginaria están separadas por un espacio, al igual que las diferentes entradas del vector resultante. Si eso parece demasiado feo, agregue un !al final ( 17 bytes ) para tener cada entrada de la salida en una fila diferente.

Explicación

-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
fuente
2

Pyth, 30

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

Banco de pruebas

Enfoque muy ingenuo, solo una implementación de la fórmula. Se encuentra con varios problemas menores de coma flotante con valores que deberían ser inversos aditivos añadidos para dar como resultado valores que están ligeramente fuera de cero.

Curiosamente .jno parece funcionar sin argumentos, pero no estoy seguro de si lo estoy usando correctamente.

FryAmTheEggman
fuente
1
¡Felicidades por 10k!
Luis Mendo
2

Pyth, 18 bytes

Utiliza el hecho de que e^ipi = -1.

La fórmula se convierte X_k = sum(x_n / ((-1)^(2nk/N))).

ms.ecb^_1*cyklQdQU

Banco de pruebas.

Monja permeable
fuente
2

Julia, 45 bytes

~=endof
!x=sum(x./im.^(4(r=0:~x-1).*r'/~x),1)

Pruébalo en línea!

El código usa la fórmula equivalente

fórmula

Dennis
fuente
2

Python 2, 78 bytes

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

El polinomio se evalúa para cada potencia pde 1j**(4./len(l)).

La expresión reduce(lambda a,b:a*p+b,l)evalúa el polinomio dado por lel valor pmediante el método de Horner:

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

Excepto, la lista de entrada se invierte, con un término constante al final. Podríamos revertirlo, pero debido p**len(l)==1a los coeficientes de Fourier, podemos usar un truco de invertir py multiplicar todo el resultado por p.

Algunos intentos de igual longitud:

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 función de 1 byte más (79):

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

Un intento de recursión (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)

Simulando iterativamente 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
fuente
2

C (gcc) , 86 78 bytes

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);}

Pruébalo en línea!

Esto supone que el vector de salida se pone a cero antes de su uso.

techo
fuente
1

Python 2, 89 bytes

Utiliza el hecho de que e^ipi = -1.

La fórmula se convierte 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!

Monja permeable
fuente
¡Publique eso como una respuesta separada!
Leaky Nun
De acuerdo, si tú lo dices.
Dennis
1

Mathematica, 49 48 47 bytes

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

Basado en la fórmula de la solución de @Dennis .

millas
fuente
1

Axioma, 81 bytes

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])

usando la fórmula que alguien publica aquí. Resultados

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