Maxima – 130 – Funzioni numeriche – 2

feynmann

Continuo da qui, copio dal Reference Manual, PDF scaricabile da qui, sono a p.384.

fft (x)
Computes the complex fast Fourier transform. x is a list or array (named or unnamed) which contains the data to transform. The number of elements must be a power of 2. The elements must be literal numbers (integers, rationals, floats, or bigfloats) or symbolic constants, or expressions a + b*%i where a and b are literal numbers or symbolic constants.

fft returns a new object of the same type as x, which is not modified. Results are always computed as floats or expressions a + b*%i where a and b are floats. If bigfloat precision is needed the function bf_fft can be used instead as a drop-in replacement of fft that is slower, but supports bfloats. In addition if it is known that the input consists of only real values (no imaginary parts), real_fft can be used which is potentially faster.

The discrete Fourier transform is defined as follows. Let y be the output of the transform. Then for k from 0 through n - 1,
y[k] = (1/n) sum(x[j] exp(+2 %i %pi j k / n), j, 0, n - 1)

As there are various sign and normalization conventions possible, this definition of the transform may differ from that used by other mathematical software.

When the data x are real, real coefficients a and b can be computed such that
x[j] = sum(a[k]*cos(2*%pi*j*k/n)+b[k]*sin(2*%pi*j*k/n), k, 0, n/2)
with
a[0] = realpart (y[0])
b[0] = 0
and, for k from 1 through n/2 - 1,
a[k] = realpart (y[k] + y[n - k])
b[k] = imagpart (y[n - k] - y[k])
and
a[n/2] = realpart (y[n/2])
b[n/2] = 0

load(fft) loads this function.

See also inverse_fft (inverse transform), recttopolar, and polartorect. See real_fft for FFTs of a real-valued input, and bf_fft and bf_real_fft for operations on bigfloat values.

(%i1) load ("fft") $

(%i2) fpprintprec : 4 $

(%i3) L : [1, 2, 3, 4, -1, -2, -3, -4] $

(%i4) L1 : fft (L);
(%o4) [0.0, 1.81 %i - 0.1035, 0.0, 0.3106 %i + 0.6035, 0.0,
                                 0.6035 - 0.3106 %i, 0.0, (- 1.81 %i) - 0.1035]
(%i5) L2 : inverse_fft (L1);
(%o5) [1.0, 4.44E-16 %i + 2.0, 3.0, 4.0 - 4.44E-16 %i, - 1.0,
                               (- 4.44E-16 %i) - 2.0, - 3.0, 4.44E-16 %i - 4.0]
(%i6) lmax (abs (L2 - L));
(%o6)                              4.44E-16

Complex data.

(%i1) load ("fft") $

(%i2) fpprintprec : 4 $

(%i3) L : [1, 1 + %i, 1 - %i, -1, -1, 1 - %i, 1 + %i, 1] $

(%i4) L1 : fft (L);
(%o4) [0.5, 0.5, 0.25 %i - 0.25, (- 0.3535 %i) - 0.3535, 0.0, 0.5,
                                        (- 0.25 %i) - 0.25, 0.3535 %i + 0.3535]
(%i5) L2 : inverse_fft (L1);
(%o5) [1.0, 1.0 %i + 1.0, 1.0 - 1.0 %i, - 1.0, - 1.0, 1.0 - 1.0 %i,
                                                             1.0 %i + 1.0, 1.0]
(%i6) lmax (abs (L2 - L));
(%o6)                                 0.0

Computation of sine and cosine coefficients.

(%i1) load ("fft") $

(%i2) fpprintprec : 4 $

(%i3) L : [1, 2, 3, 4, 5, 6, 7, 8] $

(%i4) n : length (L) $

(%i5) x : make_array (any, n) $

(%i6) fillarray (x, L) $

(%i7) y : fft (x) $

(%i8) a : make_array (any, n/2 + 1) $

(%i9) b : make_array (any, n/2 + 1) $

(%i10) a[0] : realpart (y[0]) $

(%i11) b[0] : 0 $

(%i12) for k : 1 thru n/2 - 1 do
         (a[k] : realpart (y[k] + y[n - k]),
         b[k] : imagpart (y[n - k] - y[k]));
(%o12)                               done
(%i13) a[n/2] : y[n/2] $

(%i14) b[n/2] : 0 $

(%i15) listarray (a);
(%o15)                 [4.5, - 1.0, - 1.0, - 1.0, - 0.5]
(%i16) listarray (b);
(%o16)                    [0, 2.414, 1.0, 0.4142, 0]
(%i17) f(j) := sum (a[k]*cos(2*%pi*j*k/n) + b[k]*sin(2*%pi*j*k/n), k, 0, n/2) $

(%i18) makelist (float (f (j)), j, 0, n - 1);
(%o18)             [1.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0]

real_fft (x)
Computes the fast Fourier transform of a real-valued sequence x. This is equivalent to performing fft(x), except that only the first N/2+1 results are returned, where N is the length of x. N must be power of two.

No check is made that x contains only real values.

The symmetry properites of the Fourier transform of real sequences to reduce he complexity. In particular the first and last output values of real_fft are purely real.

For larger sequencyes, real_fft may be computed more quickly than fft.

Since the output length is short, the normal inverse_fft cannot be directly used. Use inverse_real_fft to compute the inverse.

inverse_real_fft (y)
Computes the inverse Fourier transfrom of y, which must have a length of N/2+1 where N is a power of two. That is, the input x is expected to be the output of real_fft.

No check is made to ensure that the input has the correct format. (The first and last
elements must be purely real.)

bf_inverse_fft (y)
Computes the inverse complex fast Fourier transform. This is the bigfloat version of inverse_fft that converts the input to bigfloats and returns a bigfloat result.

bf_fft (y)
Computes the forward complex fast Fourier transform. This is the bigfloat version of fft that converts the input to bigfloats and returns a bigfloat result.

bf_real_fft (x)
Computes the forward fast Fourier transform of a real-valued input returning a bigfloat result. This is the bigfloat version of real_fft.

bf_inverse_real_fft (y)
Computes the inverse fast Fourier transfrom with a real-valued bigfloat output. This is the bigfloat version of inverse_real_fft.

Posta un commento o usa questo indirizzo per il trackback.

Trackback

Rispondi

Inserisci i tuoi dati qui sotto o clicca su un'icona per effettuare l'accesso:

Logo di WordPress.com

Stai commentando usando il tuo account WordPress.com. Chiudi sessione /  Modifica )

Google photo

Stai commentando usando il tuo account Google. Chiudi sessione /  Modifica )

Foto Twitter

Stai commentando usando il tuo account Twitter. Chiudi sessione /  Modifica )

Foto di Facebook

Stai commentando usando il tuo account Facebook. Chiudi sessione /  Modifica )

Connessione a %s...

Questo sito utilizza Akismet per ridurre lo spam. Scopri come vengono elaborati i dati derivati dai commenti.

%d blogger hanno fatto clic su Mi Piace per questo: