Zu schnell, zu Fourier: FFT-Code Golf

48

Implementieren Sie die schnelle Fourier-Transformation in möglichst wenigen Zeichen.

Regeln:

  • Kürzeste Lösung gewinnt
  • Es kann angenommen werden, dass der Eingang ein 1D-Array ist, dessen Länge eine Zweierpotenz ist.
  • Sie können den Algorithmus Ihrer Wahl verwenden, aber die Lösung muss tatsächlich eine schnelle Fouriertransformation sein, nicht nur eine naive diskrete Fouriertransformation (dh es müssen asymptotische Berechnungskosten von O(NLogN) ).

Bearbeiten:

  • Der Code sollte die standardmäßige schnelle Fourier-Vorwärtstransformation implementieren, deren Form in Gleichung (3) dieses Wolfram-Artikels zu sehen ist.

    Bildbeschreibung hier eingeben

  • Die Verwendung einer FFT-Funktion aus einer vorhandenen Standardbibliothek oder einem Statistikpaket ist nicht zulässig. Hier besteht die Herausforderung darin, den FFT-Algorithmus selbst kurz und bündig zu implementieren .
jakevdp
quelle
3
Dies ist unterbestimmt. Zumindest müssen Sie die Normalisierungsfaktoren definieren, und Sie sollten sich auch darüber im Klaren sein, dass Unklarheiten absichtlich falsch interpretiert werden. ZB ist "Implement" zufrieden mit der Antwort " FFT(3 Zeichen): Es ist in der Standardbibliothek"? Einige Testfälle wären auch gut.
Peter Taylor
Kommt es auf die Reihenfolge der Ausgabeelemente an, dh müssen wir bitumgekehrtes Entschlüsseln implementieren, oder können wir die Ausgabe in verschlüsselter Reihenfolge belassen?
Paul R
Siehe die Änderungen an den Regeln. Die Ausgabe sollte eine Liste / ein Array mit Werten sein, die gemäß den oben genannten Indizes im Standard-DFT-Ausdruck sortiert sind.
Jakevdp
2
Können Sie einige Beispiele für Ein- und Ausgaben veröffentlichen, damit wir unsere Implementierungen testen können?
FUZxxl
2
Titel sollte "Fast and Fourier-s" (schnell und wütend) gewesen sein.
Clismique

Antworten:

12

Mathematica, 95 Bytes

Eine weitere Implementierung der Cooley-Tukey-FFT mit Hilfe von @ chyaong .

{n=Length@#}~With~If[n>1,Join[+##,#-#2]&[#0@#[[;;;;2]],#0@#[[2;;;;2]]I^Array[-4#/n&,n/2,0]],#]&

Ungolfed

FFT[x_] := With[{N = Length[x]},
  If[N > 1,
    With[{a = FFT[ x[[1 ;; N ;; 2]] ], 
          b = FFT[ x[[2 ;; N ;; 2]] ] * Table[E^(-2*I*Pi*k/N), {k, 0, N/2 - 1}]},
      Join[a + b, a - b]],
    x]]
Meilen
quelle
1
Ich denke #[[;;;;2]]==#[[1;;N;;2]]und [[2;;;;2]]==[[2;;N;;2]].
Chyanog
1
101 Zeichen :With[{L=Length@#},If[L>1,Join[+##,#-#2]&[#0@#[[;;;;2]],#0@#[[2;;;;2]]E^(-2I*Pi(Range[L/2]-1)/L)],#]]&
Chyanog
Schön, dass Sie eine andere anonyme Funktion darin komprimieren können, ohne mit der rekursiven in Konflikt zu geraten. Außerdem wurde bekannt, dass Teil fehlende Indizes ausfüllt. Wir können es mit Unicode weiterentwickeln.
Meilen
9

J, 37 Bytes

_2&(0((+,-)]%_1^i.@#%#)&$:/@|:]\)~1<#

Eine Verbesserung nach ein paar Jahren. Verwendet weiterhin den Cooley-Tukey-FFT-Algorithmus.

4 Bytes mit e πi = -1 dank @ Leaky Nun gespeichert .

Probieren Sie es online!

Verwendungszweck

   f =: _2&(0((+,-)]%_1^i.@#%#)&$:/@|:]\)~1<#
   f 1 1 1 1
4 0 0 0
   f 1 2 3 4
10 _2j2 _2 _2j_2
   f 5.24626 3.90746 3.72335 5.74429 4.7983 8.34171 4.46785 0.760139
36.9894 _6.21186j0.355661 1.85336j_5.74474 7.10778j_1.13334 _0.517839 7.10778j1.13334 1.85336j5.74474 _6.21186j_0.355661

Erläuterung

_2&(0((+,-)]%_1^i.@#%#)&$:/@|:]\)~1<#  Input: array A
                                    #  Length
                                  1<   Greater than one?
_2&(                            )~     Execute this if true, else return A
_2                            ]\         Get non-overlapping sublists of size 2
    0                       |:           Move axis 0 to the end, equivalent to transpose
                          /@             Reduce [even-indexed, odd-indexed]
                       &$:               Call recursively on each 
                   #                     Get the length of the odd list
                i.@                      Range from 0 to that length exclusive
                    %#                   Divide each by the odd length
             _1^                         Compute (-1)^x for each x
           ]                             Get the odd list
            %                            Divide each in that by the previous
       +                                 Add the even values and modified odd values
         -                               Subtract the even values and modified odd values
        ,                                Join the two lists and return
Meilen
quelle
1
Siehe auch: blog.ndpar.com/2014/10/11/dft-j
FUZxxl
9

Python, 166 151 150 Zeichen

Hierbei wird der Cooley-Tukey-FFT-Algorithmus von Radix-2 verwendet

from math import*
def F(x):N=len(x);t=N<2or(F(x[::2]),F(x[1::2]));return N<2and x or[
a+s*b/e**(2j*pi*n/N)for s in[1,-1]for(n,a,b)in zip(range(N),*t)]

Das Ergebnis testen

>>> import numpy as np
>>> x = np.random.random(512)
>>> np.allclose(F(x), np.fft.fft(x))
True
jakevdp
quelle
1
2 Dinge: Es ist normalerweise am besten zu verwenden from x import*und sum(([x for x in y] for y in z),[])ist länger als [x for y in z for x in y].
Stand
1
Danke - das spart 15 Zeichen! 11 weitere und es ist ein Tweet.
Jakevdp
Oh, das ist definitiv möglich. Wenn Sie eine Verbesserung finden, wird eine alte oft zu einem Stolperstein.
Stand
5

Python 3: 140 134 113 Zeichen

Kurze Version - kurz und süß, passt in einen Tweet (dank Meilen ):

from math import*
def f(v):
 n=len(v)
 if n<2:return v
 a,b=f(v[::2])*2,f(v[1::2])*2;return[a[i]+b[i]/1j**(i*4/n)for i in range(n)]

(In Python 2 /wird die Division abgeschnitten, wenn beide Seiten Ganzzahlen sind. Wir ersetzen also (i*4/n)durch (i*4.0/n), wodurch die Länge auf 115 Zeichen erhöht wird.)

Lange Version - mehr Klarheit in die Interna des klassischen Cooley-Tukey FFT:

import cmath
def transform_radix2(vector):
    n = len(vector)
    if n <= 1:  # Base case
        return vector
    elif n % 2 != 0:
        raise ValueError("Length is not a power of 2")
    else:
        k = n // 2
        even = transform_radix2(vector[0 : : 2])
        odd  = transform_radix2(vector[1 : : 2])
        return [even[i % k] + odd[i % k] * cmath.exp(i * -2j * cmath.pi / n) for i in range(n)]
Nayuki
quelle
1
Verkürzt auf 113 Bytes mite^(-2j * pi * i / n) = (-1)^(2 * i / n) = (1j)^(4 * i / n)
Meilen
@miles Sehr beeindruckende Beobachtung, danke! Nachdem ich DFTs über ein Jahrzehnt lang wiederholt implementiert hatte, war ich von sin / cos / exp besessen und vergaß, dass einfache Potenzen von i verwendet werden können. Ich habe meine Antwort bearbeitet, um die neuen Erkenntnisse einzubeziehen und Sie zu würdigen.
Nayuki
5

R: 142 133 99 95 Bytes

Vielen Dank an @ Giuseppe für die Hilfe beim Rasieren von 32 bis 36 Bytes!

f=function(x,n=sum(x|1),y=1:(n/2)*2)`if`(n>1,f(x[-y])+c(b<-f(x[y]),-b)*exp(-2i*(y/2-1)*pi/n),x)

Ein zusätzlicher Trick besteht darin, die Standardargumente der Hauptfunktion zu verwenden, um einige Variablen zu instanziieren.
Die Verwendung ist immer noch die gleiche:

x = c(1,1,1,1)
f(x)
[1] 4+0i 0+0i 0+0i 0+0i

4 Jahre alte Version mit 133 Bytes:

f=function(x){n=length(x);if(n>1){a=Recall(x[seq(1,n,2)]);b=Recall(x[seq(2,n,2)]);t=exp(-2i*(1:(n/2)-1)*pi/n);c(a+b*t,a-b*t)}else{x}}

Mit Einrückungen:

f=function(x){
    n=length(x)
    if(n>1){
        a=Recall(x[seq(1,n,2)])
        b=Recall(x[seq(2,n,2)])
        t=exp(-2i*(1:(n/2)-1)*pi/n)
        c(a+b*t,a-b*t)
        }else{x}
    }

Es wird auch der Cooley-Tukey-Algorithmus verwendet. Die einzigen Tricks hier sind die Verwendung von Funktionen Recall, die Rekursivität ermöglichen, und die Verwendung von R-Vektorisierung, die die tatsächliche Berechnung erheblich verkürzen.

Verwendungszweck:

x = c(1,1,1,1)
f(x)
[1] 4+0i 0+0i 0+0i 0+0i
Plannapus
quelle
1
Vier Jahre später erreichen wir 101 Bytes . Nicht 100% ig sicher, warum Sie Recalldiese Funktion verwendet haben, aber hey, es ist im Nachhinein einfach, Golf zu spielen! :) +1, sehr nett.
Giuseppe
Ja Recallist jetzt in der Tat unnötig. Ich habe das vor ein paar Monaten bemerkt, war aber zu faul, um es zu ändern :) Ich werde es ändern.
Plannapus
Sehr schön! Ich habe weitere 4 Bytes herausgepresst! und stellt dies Mathematica gleich.
Giuseppe
Vielen Dank! Ich habe darüber nachgedacht, mich dort niederzulassen y, aber ich habe nicht bemerkt, dass es auch für das exp(...)Teil verwendet werden könnte.
Plannapus
4

Python, 134

Dies basiert stark auf der Lösung von jakevdp, daher habe ich dieses Wiki auf ein Community-Wiki gesetzt.

from math import*
F=lambda x:x*(len(x)<2)or[a+s*b/e**(2j*pi*n/len(x))for s in(1,-1)for n,(a,b)in
enumerate(zip(F(x[::2]),F(x[1::2])))]

Änderungen:

-12 Zeichen: töten t.

def F(x):N=len(x);t=N<2or(F(x[::2]),F(x[1::2]));return ... in zip(range(N),*t)]
def F(x):N=len(x);return ... in zip(range(N),F(x[::2]),F(x[1::2]))]

-1 Zeichen: Exponententrick x*y**-z == x/y**z (dies könnte einigen anderen helfen)

...[a+s*b*e**(-2j*pi*n/N)...
...[a+s*b/e**(2j*pi*n/N)...

-2 Zeichen: Ersetzen anddurch*

...return N<2and x or[
...return x*(N<2)or[

+1 lambdaZeichen : Ize, tötenN

def F(x):N=len(x);return x*(N<2)or[a+s*b/e**(2j*pi*n/N) ... zip(range(N) ...
F=lambda x:x*(len(x)<2)or[a+s*b/e**(2j*pi*n/len(x)) ... zip(range(len(x)) ...

-2 Zeichen: Verwenden Sie enumerateanstelle vonzip(range(len(

...for(n,a,b)in zip(range(len(x)),F(x[::2]),F(x[1::2]))]
...for n,(a,b)in enumerate(zip(F(x[::2]),F(x[1::2])))]
boothby
quelle
Ich denke, dass dies keine schnelle Fouriertransformation mehr ist. Indem Sie t "töten", haben Sie einige unnötige Berechnungen hinzugefügt, die es von O [N log (N)] nach O [N ^ 2] verschieben
jakevdp
Es scheint, dass ich meinen eigenen Beitrag nicht ablehnen kann. Sie haben Recht, ich habe die Loop-Reihenfolge vertauscht und die Aufführung abgebrochen. Ich lasse dies für den Fall, dass ich einen Weg finde, es zu beheben.
Stand
101 Bytes mitf=lambda x:x*(len(x)<2)or[u+v/1j**(4*i/len(x))for i,(u,v)in enumerate(zip(f(x[::2])*2,f(x[1::2])*2))]
Meilen
Sie können ersetzen for s in(1,-1)formit for s in 1,-1foroder sogar, wenn Reihenfolge keine Rolle spielt, for s in-1,1for.
Jonathan Frech
4

C 259

typedef double complex cplx;
void fft(cplx buf[],cplx out[],int n,int step){
if(step < n){
fft(out, buf,n, step * 2);
fft(out+step,buf+step,n,step*2);
for(int i=0;i<n;i+=2*step){
cplx t=cexp(-I*M_PI*i/n)*out[i+step];
buf[i/2]=out[i]+t;
buf[(i+n)/2]=out[i]-t;
}}}

Das Problem ist, dass solche Implementierungen nutzlos sind und der einfache Algorithmus VIEL schneller ist.

Sanaris
quelle
2
Sie können mehr Leerzeichen entfernen, um eine geringere Anzahl von Zeichen zu erhalten. Beispielsweise step < nkann in step<nund step * 2in geändert werden step*2.
ProgramFOX
2
Alle Variablen, Funktionen und typedefs sollten Namen mit einem Buchstaben haben, um viele Zeichen zu speichern
2
Sie haben von jemandem eine Reihe von Verbesserungen vorschlagen lassen. Schau sie dir
Justin
1
Sie können alle Zeilenumbrüche entfernen, Zeilenumbrüche sind in C
TuxCrafting am
@ TùxCräftîñg Nicht alle Zeilenumbrüche sind unbrauchbar. Sie werden für Direktiven wie #include, #define, #if usw. benötigt
Nayuki
3

Matlab, 128 118 107 102 101 94 93 Bytes

EDIT6: danke @algmyr für ein weiteres Byte!

function Y=f(Y);
n=numel(Y);
k=2:2:n;
if k;
   c=f(Y(k-1));
   d=f(Y(k)).*i.^(2*(2-k)/n);
   Y=[c+d;c-d];
end

EDIT5: Immer noch kürzer :) dank @sanchises

function Y=f(Y)
n=numel(Y);
k=2:2:n;
if k;
   c=f(Y(k-1));
   d=f(Y(k)).*(-1).^((2-k)/n);
   Y=[c+d;c-d];
end

EDIT4: Yay, -1 Zeichen mehr (hätte auch ohne das gehen können k):

function Y=f(Y)
n=numel(Y);
if n>1;
   k=2:2:n;
   c=f(Y(k-1));
   d=f(Y(k)).*(-1).^((k/2-1)*2/n)';
   Y=[c+d;c-d];
end

EDIT2 / 3: Danke für @sanchises für weitere Verbesserungen!

function Y=f(Y)
n=numel(Y);  
if n>1;
   c=f(Y(1:2:n));
   d=f(Y(2:2:n)).*(-1).^(-(0:n/2-1)*2/n).';
   Y=[c+d;c-d]; 
end

BEARBEITEN: Konnte einige Verbesserungen vornehmen und stellte fest, dass die Skalierungskonstante nicht erforderlich ist.

Dies ist die erweiterte Version. Die Anzahl der Zeichen ist gültig, wenn Sie die Zeilenumbrüche / Leerzeichen entfernen. (Funktioniert nur für Spaltenvektoren.)

function y=f(Y)
n=numel(Y);  
y=Y;
if n>1;
   c=f(Y(1:2:n));
   d=f(Y(2:2:n));
   n=n/2;
   d=d.*exp(-pi*i*(0:n-1)/n).';
   y=[c+d;c-d]; 
end
fehlerhaft
quelle
Tipp: Sie können die beiden kombinieren d=Linien in einem: m=n/2;d=f(Y(2:2:n)).*exp(-pi*i*(0:m-1)/m).';. Darüber hinaus halten Wechsel y=f(Y)zu Y=f(Y)und entfernen Linie 3 (und versprechen Sie nie , dass außerhalb von Code-Golf tun)
Sanchises
Oh Danke! Hat function Y = f(Y)es andere Nachteile als Unlesbarkeit?
Fehler
Nun, MATLAB wird sich nie über einen Rückgabewert beschweren, auch wenn Y nie geändert wird. Es ist zwar etwas schneller, also denke ich, dass es für einige Zwecke doch nicht so schlimm ist (dh eine Funktion, die die Eingabevariable fast nie ändert)
Sanchises
Jetzt, um mehr zu rasieren: m=n/2könnte entfernt und stattdessen mdurch n/2und n*2ersetzt werden. Und dann glaube ich fest daran, dass das Programm so kurz ist, wie es in MATLAB sein könnte.
Sanchises
1
Und dann glaube ich fest daran, dass das Programm so kurz ist, wie es in MATLAB sein könnte. - Sanchises 8. März 15 um 21:05 Uhr Berühmte letzte Worte ...
Sanchises
2

Jelly, 31 30 28 26 Bytes , nicht konkurrierend

LḶ÷$N-*×,N$+ḷF
s2Z߀ç/µ¹Ṗ?

Jelly wurde nach dieser Herausforderung entwickelt und ist daher nicht konkurrierend.

Dies verwendet den rekursiven Cooley-Tukey-Radix-2-Algorithmus. Für eine Version ohne Golf siehe meine Antwort in Mathematica.

Probieren Sie es online aus oder überprüfen Sie mehrere Testfälle .

Erläuterung

LḶ÷$N-*×,N$+ḷF  Helper link. Input: lists A and B
L               Get the length of A
   $            Operate on that length
 Ḷ                Make a range [0, 1, ..., length-1]
  ÷               Divide each by length
    N           Negate each
     -          The constant -1
      *         Compute -1^(x) for each x in that range
       ×        Multiply elementwise between that range and B, call it B'  
          $     Operate on that B'
         N        Negate each
        ,         Make a list [B', -B']
            ḷ   Get A
           +    Add vectorized, [B', -B'] + A = [A+B', A-B']
             F  Flatten that and return

s2Z߀ç/µ¹Ṗ?  Main link. Input: list X
         Ṗ   Curtail - Make a copy of X with the last value removed
          ?  If that list is truthy (empty lists are falsey)
       µ       Parse to the left as a monad
s2             Split X into sublists of length 2
  Z            Transpose them to get [even-index, odd-index]
   ߀          Call the main link recursively on each sublist
     ç/        Call the helper link as a dyad on the sublists and return
             Else
        ¹      Identity function on X and return
Meilen
quelle
2

C (gcc) , 188 186 184 183 Bytes

#define d(a,b,c)f(a,b,c,1,0)
f(a,b,c,n,k)_Complex*a,*b;{_Complex z[c];*b=*a;if(n<c)for(f(a,z,c,n*2),f(a+n,z+n,c,n*2);k<c;k+=n*2)b[k+c>>1]=z[k]*2-(b[k/2]=z[k]+z[k+n]/cpow(1i,2.*k/c));}

Probieren Sie es online!

Etwas weniger golfen

#define d(a,b,c)f(a,b,c,1,0)
f(a,b,c,n,k)_Complex*a,*b;{
  _Complex z[c];
  *b=*a;
  if(n<c)
    for(f(a,z,c,n*2),f(a+n,z+n,c,n*2);k<c;k+=n*2)
      b[k+c>>1]=z[k]*2-(b[k/2]=z[k]+z[k+n]/cpow(1i,2.*k/c));
}
Ceilingcat
quelle
1

Pari / GP, 76 Zeichen

X(v)=my(t=-2*Pi*I/#v,s);vector(#v,k,s=t*(k-1);sum(n=0,#v-1,v[n+1]*exp(s*n)))

Verwendungszweck

X([1,1,1,1])
%2 = [4.000000000000000000000000000, 0.E-27 + 0.E-28*I, 0.E-28 + 0.E-27*I, 0.E-27 + 0.E-28*I]
P̲̳x͓L̳
quelle
3
Ist das nicht die naive DFT? (dh Theta (N ^ 2))
Meilen
1

Oktave , 109 103 101 100 Bytes

f(f=@(f)@(x,n=rows(x)){@(d=f(f)(x(k=2:2:n)).*i.^((k*2-4)/n)')[d+(c=f(f)(x(k-1)));c-d],x}{1+(n<2)}())

Probieren Sie es online!

Ooooo bluten meine Augen von diesem rekursiven verfluchten Lambda. Große Teile davon wurden aus der Antwort von @flawr gestrichen.

f(                                          % lambda function
  f=@(f)                                    % defined in its own argument list, 
                                            % accepts itself as parameter (for recursion)
    @(x,n=rows(x)){                         % calls another lambda,
                                            % 2nd parameter is just to define a variable
      @(d=f(f)(x(k=2:2:n)).*i.^((k*2-4)/n)')% 1/4 of FFT (argument just defines a variable)
        [d+(c=f(f)(x(k-1)));                % 2/4 of FFT
         c-d                                % 4/4 of FFT
        ],                                  % This is in a @()[] to inhibit evaluation
                                            % unless actually called
      x                                     % FFT of length 1
    }{1+(n<2)}                              % if len(x)==1, return x
                                            % else return [d+c;c-d]
  ()                                        % this is probably important too
)
Ceilingcat
quelle
Ich verstehe nicht, was du getan hast, aber ich mag es sehr.
Fehler
0

Axiom, 259 , 193 , 181 , 179 Bytes

L(g,n,f)==>[g for i in 1..n|f]
h(a)==(n:=#a;n=1=>a;c:=h(L(a.i,n,odd? i));d:=h(L(a.i,n,even? i));n:=n/2;t:=1>0;v:=L(d.i*%i^(-2*(i-1)/n),n,t);append(L(c.i+v.i,n,t),L(c.i-v.i,n,t)))

Selbst wenn h (a) den gesamten Test bestehen könnte und als Eintrag für diesen "Wettbewerb" in Ordnung wäre, muss man unten h () oder hlp () bis fft () aufrufen, um die Argumente zu überprüfen . Ich weiß nicht, ob diese Software in Ordnung ist, weil ich nur gesehen habe, was andere geschrieben haben, und nach der Art und Weise gesucht habe, wie sie in Axiom ausgeführt werden kann, um ein mögliches richtiges Ergebnis zu erhalten. Unter ungolfed Code mit ein paar Kommentaren:

-- L(g,n,f)==>[g for i in 1..n|f]
-- this macro L, build one List from other list, where in g, there is the generic element of index i
-- (as a.i, or a.i*b.i or a.i*4), n build 1..n that is the range of i, f is the condition 
-- for insert the element in the list result.

hlp(a)==
    n:=#a;n=1=>a
    -- L(a.i,n,odd? i)  it means build a list getting "even indices i of a.i as starting from index 0" [so even is odd and odd is even]
    -- L(a.i,n,even? i) it means build a list getting "odd  indices i of a.i as starting from index 0"
    c:=hlp(L(a.i,n,odd? i));d:=hlp(L(a.i,n,even? i))
    n:=n/2;t:=1>0
    v:=L(d.i*%i^(-2*(i-1)/n),n,t)
    append(L(c.i+v.i,n,t),L(c.i-v.i,n,t))

-- Return Fast Fourier transform of list a, in the case #a=2^n
fft(a)==(n:=#a;n=0 or gcd(n,2^30)~=n=>[];hlp(a))

(5) -> h([1,1,1,1])
   (5)  [4,0,0,0]
                                    Type: List Expression Complex Integer
(6) -> h([1,2,3,4])
   (6)  [10,- 2 + 2%i,- 2,- 2 - 2%i]
                                    Type: List Expression Complex Integer
(7) -> h([5.24626,3.90746,3.72335,5.74429,4.7983,8.34171,4.46785,0.760139])
   (7)
   [36.989359, - 6.2118552150 341603904 + 0.3556612739 187363298 %i,
    1.85336 - 5.744741 %i, 7.1077752150 341603904 - 1.1333387260 812636702 %i,
    - 0.517839, 7.1077752150 341603904 + 1.1333387260 812636702 %i,
    1.85336 + 5.744741 %i,
    - 6.2118552150 341603904 - 0.3556612739 187363298 %i]
                                      Type: List Expression Complex Float
(8) -> h([%i+1,2,%i-2,9])
   (8)  [10 + 2%i,3 + 7%i,- 12 + 2%i,3 - 7%i]
                                    Type: List Expression Complex Integer

In den wenigen Fällen, in denen ich h () oder fft () gesehen hatte, würde die exakte Lösung zurückgegeben, aber wenn die Vereinfachung nicht gut ist, wie in:

(13) -> h([1,2,3,4,5,6,7,8])
   (13)
                    +--+                                   +--+
        (- 4 + 4%i)\|%i  - 4 + 4%i             (- 4 - 4%i)\|%i  - 4 + 4%i
   [36, --------------------------, - 4 + 4%i, --------------------------, - 4,
                    +--+                                   +--+
                   \|%i                                   \|%i
            +--+                                   +--+
    (- 4 + 4%i)\|%i  + 4 - 4%i             (- 4 - 4%i)\|%i  + 4 - 4%i
    --------------------------, - 4 - 4%i, --------------------------]
                +--+                                   +--+
               \|%i                                   \|%i
                                    Type: List Expression Complex Integer

als es genug ist, ändern Sie den Typ von nur einem Element der Liste, wie im folgenden Schreiben von 8. (Float), um die ungefähre Lösung zu finden:

(14) -> h([1,2,3,4,5,6,7,8.])
   (14)
   [36.0, - 4.0000000000 000000001 + 9.6568542494 923801953 %i, - 4.0 + 4.0 %i,
    - 4.0 + 1.6568542494 92380195 %i, - 4.0, - 4.0 - 1.6568542494 92380195 %i,
    - 4.0 - 4.0 %i, - 4.0 - 9.6568542494 923801953 %i]
                                      Type: List Expression Complex Float

Ich habe es geschrieben, alle anderen Antworten gesehen, weil es in dem Link auf der Seite zu schwierig war, so dass ich nicht weiß, ob dieser Code richtig sein kann. Ich bin kein FFT-Experte, also kann all dies (es ist wahrscheinlich) falsch sein.

RosLuP
quelle
0

APL (NARS), 58 Zeichen, 116 Byte

{1≥k←≢⍵:⍵⋄(∇⍵[y∼⍨⍳k])(+,-)(∇⍵[y←2×⍳t])×0J1*t÷⍨2-2×⍳t←⌊k÷2}

Prüfung

  f←{1≥k←≢⍵:⍵⋄(∇⍵[y∼⍨⍳k])(+,-)(∇⍵[y←2×⍳t])×0J1*t÷⍨2-2×⍳t←⌊k÷2}
  f 1 1 1 1
4J0 0J0 0J0 0J0 
  f 1 2 3 4
10J0 ¯2J2 ¯2J0 ¯2J¯2 
  f 1J1 2 ¯2J1  9
10J2 3J7 ¯12J2 3J¯7 
  f 5.24626,3.90746,3.72335,5.74429,4.7983,8.34171,4.46785,0.760139
36.989359J0 ¯6.211855215J0.3556612739 1.85336J¯5.744741 7.107775215J¯1.133338726 ¯0.517839J0 
  7.107775215J1.133338726 1.85336J5.744741 ¯6.211855215J¯0.3556612739 
RosLuP
quelle