ares@alessia.dei.unipd.it (Nicola Catacchio 259126) (05/09/91)
I wrote a routine for fft.Accept one argument,a real or complex
vector whose size must be a power of 2.
It's about 33% faster than a 'very fft' program that was on the news
some time ago.To make a representation of the power spectrum of a function,
use SAMP: it takes a function from level 2 and a list containing two real
numbers that define the sampling interval, and a variable appearing in
the function in level 2. The global variable 'n', already defined in the
directory, must contain the # of samples, a power of 2 ( 64 ).
The procedure SAMP produces a real vector that can be transformed
using FFT, to get the discrete Fourier transform.
The complex vector left on the stack by FFT being complex, cannot
be plotted as it is: Using ABSV, it will be reduced at his absolute value,
that is, his power component.
Now use PLOT : a histogram will be plotted on the screen, and to
save memory, when you exit from GRAPH, it erases the PICT and PPAR.
Example: take the power spectrum of the fuction 'SIN(X)', sampled
between -20 and 20.
'SIN(X)' [ENTER]
{-20 20 X} [ENTER] [SAMP] [FFT] [ABSV] [PLOT]
You'll get a histogram representing two peaks very close one to
the other.
Sorry for any grammar mistakes: pleae E-mail to me any error
or malfunction you eventually find in programs.
|------------------------------------------------------------------------------|
|Nicola Catacchio |E-mail: ares@alessia.unipd.it |
|Universita' di Padova |mail : Cannaregio 4389, Venezia, Italy |
| |phone#: 041/5222516 |
|------------------------------------------------------------------------------|
------------------------CUT HERE------------------------------------------------
%%HP: T(3)A(R)F(.);
DIR
REVRS
\<< OBJ\-> EVAL \-> N
\<<
IF N 2 >
THEN 2 N 1
-
FOR I I
ROLL
NEXT
END N \->ARRY
\>>
\>>
FFFT
\<< OBJ\-> EVAL
DUBL IBIT DTEM
\->ARRY 2 / DUP REVRS
DUP2 - ROT ROT +
C\->R ROT C\->R SWAP
NEG SWAP ROT SWAP
R\->C CPT ROT ROT R\->C
DUP2 - ROT ROT +
\>>
CPT
\<< OBJ\-> EVAL DUP
1 + \-> N M
\<< \pi i * \->NUM
N / NEG EXP 1 3 N 1
+
FOR I OVER
* I ROLL OVER * I
ROLLD
NEXT DROP2
N \->ARRY
\>>
\>>
DUBL
\<< \-> N
\<< N DUP 2 / 1
+
FOR I I
ROLL I ROLL SWAP
R\->C -1
STEP N 2 /
\>>
\>>
SAMP
\<< EVAL 3 PICK 3
PICK SWAP - n / \-> F
A B X S
\<< A X STO 1 n
START F
EVAL S X STO+
NEXT n
\->ARRY
\>>
\>>
n 64
ABSV
\<< OBJ\-> EVAL \-> N
\<< 1 N
START N
ROLL C\->R SQ SWAP SQ
+
NEXT N
\->ARRY
\>>
\>>
PLOT
\<< DUP OBJ\-> EVAL
\-> N
\<< 1 N 2 /
START N
ROLL
NEXT { N 1
} \->ARRY STO\GS
BARPLOT GRAPH {
PPAR \GSPAR \GSDAT PICT
} PURGE
\>>
\>>
IFFT
\<< OBJ\-> EVAL \-> N
\<< 2 N 1 -
FOR I I
ROLL
NEXT N IBIT
DTEM \->ARRY N /
\>>
\>>
FFT
\<< OBJ\-> 1 GET
IBIT DTEM \->ARRY
\>>
DTEM
\<< 1 OVER \-> H N
\<< -1 1 ROT LN
2 LN /
START 1 N 2
+ DUP H - 1 +
FOR J J 3
FOR I I
H - ROLL OVER * I
ROLL SWAP DUP2 -
ROT ROT + I ROLLD I
H - ROLLD H DUP +
NEG
STEP
OVER * -1
STEP DROP
H DUP + 'H' STO
CONJ \v/ CONJ
NEXT DROP N
\>>
\>>
IBIT
\<<
IF DUP 2 >
THEN DUP 2 /
DUP \-> N
\<< 1 + 3 ROT
FOR I
IF DUP
I 1 - >
THEN I
ROLL OVER 1 + ROLLD
DUP ROLL I ROLLD
END N
SWAP
WHILE
DUP2 <
REPEAT
OVER - SWAP 2 /
SWAP
END +
NEXT
\>>
END
\>>
END