mcgrant@elaine3.stanford.edu (Michael Grant) (11/22/90)
Many of you asked for my programs concerning 1) matrix resolvents,
and 2) DFT's and convolutions of discrete functions. The former
shall come at a future time, but, for now, here are the DFT
programs.
All of these functions operate on column vectors. Originally,
they worked on lists, but there are so many functions that
can be performed on the elements of an array that I just HAD
to convert them. SO, they don't follow conventional notation
like they used to, but they are now part of a much more
powerful repertoire of functions.
The following programs are a set that I have compiled for my
OWN benefit to perform some simple DFT analysis for academic
functions. They are by no means awesome; they only handle one
dimension, and they aren't terribly fast. BUT, they sure have
saved time on homework and exams...
While I THINK they are bug-free, I had to type them in by hand,
so who knows. Of course, I don't musspelll tu muche when I type,
so perhaps these will be OK. In any case, I would appreciate it
if you notify me of any corrections that need to be made
(at mcgrant@portia.stanford.edu).
NOTE: I will post an un-commented set of these next, so that
you will have less editing to do. Note that I have a 28S, so
these programs may have to be tweaked...
---
I recommend this order in your directory:
{ DFT IDFT DHT IDHT CCNV LCNV PROD MP->C C->MP
PAD2 RPAD ZPAD SAMP AR01 AR02 PRCSN } ORDER
'DFT' -- performs a Discrete Fourier Transform
on the given input vector.
INPUT: 1: a column vector L
OUTPUT: 1: DFT(L)
REQUIRES: DHT
COMMENTS: This program uses the Discrete Hartley Transform.
The DHT is a lot faster than the DFT, and the
time required to extract the DFT from the DHT is
minimal in comparison. Of course, for real speed,
perhaps someday I will write an FHT/FFT!
<< DHT DUP ARRY-> SZ
<< 2 SZ LIST-> DROP 1 -
FOR I I ROLL
NEXT SZ ->ARRY
>> DUP2 - (0,1) * - +
IF DUP IM CNRM NOT
THEN RE
END .5 *
>>
'IDFT' -- performs an Inverse Discrete Fourier Transform
INPUT: 1: a column vector L
OUTPUT: 1: IDFT(L)
REQUIRES: DFT
COMMENTS: Nothing more than INV(N)*CONJ(DFT(CONJ(L)).
<< CONJ DFT CONJ DUP SIZE LIST-> DROP / >>
'DHT' -- performs a Discrete Hartley Transform
INPUT: 1: a column vector L
OUTPUT: 1: DHT(L)
REQUIRES: PRCSN
COMMENTS: The Discrete Hartley Transform uses cas(theta)=
cos(theta)+sin(theta) as its kernel. Some people
will notice that this can be calculated as sqrt(2)*
sin(theta+pi/4), but for some reason this calculation
just wasn't quite as accurate as sin+cos; besides,
the speed savings is not all that great.
The PRCSN variable tells the DHT to round the answers
to PRCSN digits. I do this because 4.0000000001 takes
up the whole screen, tho' I know it's 4!
Notice that ~pi~ is the pi key .
<< DUP SIZE LIST-> DROP RCLF -> L N F
<< RAD PRCSN FIX L ~pi~ ~pi~ + -> S
<< 0 N 1 -
FOR J J N / S * DUP SIN SWAP COS + J 1 + SWAP PUT
NEXT
>> -> C
<< L 0 N 1 -
FOR J 0 0 N 1 -
FOR K J * N MOD 1 + GET 'L' K 1 + GET * +
NEXT RND J 1 + SWAP PUT
NEXT
>> F STOF
>>
>>
'IDHT' -- performs an Inverse Discrete Hartley Transform.
INPUT: 1: a column vector L
OUTPUT: 1: IDHT(L)
REQUIRES: DHT
COMMENTS: DHT is symmetric, so this is easy (N^-1*DHT(L)).
<< DHT DUP SIZE LIST-> DROP / >>
'SAMP' -- generates an array by sampling a function.
INPUT: 2: a self-contained function F
1: the number of samples to take
OUTPUT: [ F(0) F(1) F(2) ... F(N-1) ]
COMMENTS: Sorry this isn't more complex guys and gals, but
it works fine for me.
<< -> F N
<< 0 N 1 -
FOR J J F EVAL
NEXT { N } ->ARRY
>>
>>
'ARO2' -- performs a binary operation on the elements of
two vectors (thereby combining them).
INPUT: 3: column vector L1 (size S1)
2: column vector L2 (size S2)
1: a self-contained binary function F
OUTPUT: [ F(L1(1),L2(1)) F(L1(2),L2(2)) ...
F(L1(min(S1,S2)),L2(min(S1,S2))) ]
COMMENTS: Notice that, if the two lists are of different
sizes, that AR02 truncates (in effect) the
longer one to match the size of the shorter one.
<< -> L1 L2 PR
<< L1 SIZE LIST-> DROP L2 SIZE LIST-> DROP MIN -> N
<< L1 { N } RDM 1 N
FOR J 'L1' J GET 'L2' J GET PR EVAL J SWAP PUT
NEXT
>>
>>
>>
'PROD' -- Performs an element-by-element multiplication.
INPUT: 2: column vector L1 (size S1)
1: column vector L2 (size S2)
OUTPUT: [ L1(1)*L2(1) L1(2)*L2(2) ... L1(min(S1,S2))*L2(min(S1,S2)) ]
REQUIRES: ARO2
COMMENTS: See ARO2. This function is useful because convolution
in one domain is equivalent to multiplication in the
other domain...but this is not QUITE true for Hartley...
<< << * >> ARO2 >>
'AR01' -- Performs a function on each element of the vector.
INPUT: 2: column vector L (size N)
1: self-contained function F
OUTPUT: 1: [ F(1) F(2) ... F(N) ]
COMMENTS: well, I just use this for C->MP and MP->C below.
<< -> L PR
<< L SIZE LIST-> DROP -> N
<< L 1 N
FOR J 'L' J GET PR EVAL J SWAP PUT
NEXT
>>
>>
>>
'C->MP' -- converts a vector from complex to mag/phase.
INPUT: 1: column vector L (size N)
OUTPUT: 1: [ R->P(L(1)) R->P(L(2)) ... R->P(L(N)) ]
REQUIRES: AR01
COMMENTS: Convenient when you just want the magnitude
of your DFT (just do a C->MP RE)!
<< (1,0) * << R->P >> ARO1 >>
'MP->C' -- converts a vector from mag/phase to complex.
INPUT: 1: column vector L (size N)
OUTPUT: 1: [ P->R(L(1)) P->R(L(2)) ... P->R(L(N)) ]
REQUIRES: AR01
COMMENTS: When you're too short-sighted to remember to
save the original answer to the DFT...
<< (1,0) * << P->R >> ARO1 >>
'CCNV' -- performs the circular convolution of two vectors
INPUT: 2: column vector L1 (size N)
1: column vector L2 (size N)
OUTPUT: 1: L1 C* L2 (size N)
COMMENTS: This program will generate an error if the two lists
are of different sizes. So, use ZPAD, PAD2, or
RPAD to adjust the sizes of your vectors.
<< DUP SIZE LIST-> DROP -> L1 L2 N
<< IF N L1 SIZE LIST-> DROP <>
THEN [ 1 ] TRN
ELSE L1 1 N
FOR J 0 1 N
FOR K 'L1' K GET 'L2' J K - N MOD 1 + GET * +
NEXT J SWAP PUT
NEXT
END
>>
>>
'LCNV' -- performs the linear convolution of two vectors
INPUT: 2: column vector L1 (size N1)
1: column vector L2 (size N2)
OUTPUT: 1: L1 * L2 (size N1+N2-1)
REQUIRES: PAD2 CCNV
COMMENTS: Gee, the REQUIRES line is longer than the program!
<< PAD2 CCNV >>
'PAD2' -- pads two vectors so L1 C* L2 <==> L1 * L2
INPUT: 2: column vector L1 (size N1)
1: column vector L2 (size N2)
OUTPUT: 2: column vector L1 (size N1+N2-1)
1: column vector L2 (size N1+N2-1)
COMMENTS: This program is essential to the proper operation of
LCNV, because circular convolution must 'see' enough
padded zeros so that no circular overlap occurs.
<< DUP2 SIZE LIST-> DROP SWAP SIZE LIST-> DROP + 1 - 1 ->LIST -> S
<< S RDM SWAP S RDM SWAP >>
>>
'ZPAD' -- adds as many zeros as is necessary to extend the
input vector to at least N elements
INPUT: 2: column vector L (size M)
1: the desired new vector size N
OUTPUT: IF N>M, 1: RDM(L,N)
IF N<M, 2: L
COMMENTS: Not much to this.
<< -> L N
<< L
IF L SIZE LIST-> DROP N <
THEN { N } RDM
END
>>
>>
'RPAD' -- replicate the vector to hit N vectors.
INPUT: 2: column vector L (size M)
1: the desired new vector size N
OUTPUT: IF N>M: [ L(1) L(2) ... L(M) L(1) .. L(M) ... ] (size N)
IF N<M: L
COMMENTS: This is useful when circularly convolving two functions
one's period is a multiple of the other's. Otherwise,
I can't think of any mathematically safe use for this.
<< -> L N
<< L SIZE LIST-> DROP -> M
<< L
IF M N <
THEN { N } RDM M 1 + N
FOR J 'L' J M MOD
IF DUP NOT
THEN DROP M
END GET J SWAP PUT
NEXT
END
>>
>>
>>
'PRCSN' -- # of digits to RND to in DHT
COMMENTS: I recommend a fairly large number, though my
experience is that the DHT is accurate to about
10 digits (as far as roundoff error goes).
10
---
Michael C. Grant "God does not play dice." Einstein
Information Systems Lab "Geez, He'd win a lot if he did,
mcgrant@portia.stanford.edu though..." Mike