ryoder@ecn.purdue.edu (Robert W Yoder) (07/28/90)
Here is a set of programs for digital signal processing.
Some have been posted before, and some are modifications of prior postings.
The DFT and FFT are the fastest/smallest I've found.
The Two-Dimensional Fourier programs are all new.
These are posted as a set because many of them call others in the set.
The following programs are called, but not listed here:
"FAST" is a program thats sets the HP28 to maximum speed.
"ROUN" is a program that cleans up the output to eliminate the inherent
computational errors involving fractional parts containing trailing
9's, or infinitesimal values. I highly recommend using it.
Included are:
Discrete Fourier Transform (also does inverse)
Fast Fourier Transform
Inverse Fast Fourier Transform
Fourier Transform (selects FFT or DFT depending on size)
Inverse Fourier Transform (selects FFT or DFT depending on size)
Periodic Convolution
Aperiodic Convolution
Two Dimensional Fourier Transform
Inverse Two Dimensional Fourier Transform
Axis Shifter for use with Two Dimensional Fourier Transforms
--------------------------------------------------------------------------------
Discrete Fourier Transform
input: level 2; vector
input: level 1; '1' for for fourier transform, '-1' for inverse transform
output:level 1; transformed vector
DFT [825E]
<< SWAP
DUP SIZE 1 GET [pi] ->NUM (0,-2) * OVER / 4 PICK * -> d x s q
<< 0 s 1 - FOR k 0
0 s 1 - FOR n q k * n * EXP x n 1 + GET * + NEXT
NEXT
s ->ARRY
IF d -1 SAME
THEN s INV *
END
>> >>
--------------------------------------------------------------------------------
Fast Fourier Transform by lehmann_a%ELDE.EPFL.CH@VM1.NoDak.EDU
input: level 2; vector
output:level 1; transformed vector
FFT [F582]
<< ARRY-> 1 GET IF DUP LN 2 LN / DUP IP == THEN
-> t
<< 1 t LN 2 LN / FOR c -2 [pi] * i * 2 c ^ / EXP ->NUM
2 c 1 - ^ -> f k
<< t 2 / 1 + t FOR d 0 k 1 - FOR e t ROLL d e + ROLL f e
^ * - LAST + e 2 + ROLLD NEXT k STEP
>>
NEXT t ->ARRY
>>
ELSE ->ARRY "FFT Error:
Size <> 2^N" 1 DISP 1400 .07 BEEP
END
>>
--------------------------------------------------------------------------------
Inverse Fast Fourier Transform by: lehmann_a%ELDE.EPFL.CH@VM1.NoDak.EDU
input: level 2; vector
output:level 1; transformed vector
FFTI [F7EE]
<< CONJ FFT DUP SIZE 1 GET / CONJ >>
--------------------------------------------------------------------------------
Fourier Transform
FT checks size and calls FFT if size is a power of 2, else calls DFT.
input: level 1; vector
output:level 1; transformed vector
FT [D998]
<< FAST DUP SIZE 1 GET
IF LN 2 LN / FP
THEN 1 DFT
ELSE FFT
END
ROUN
>>
--------------------------------------------------------------------------------
Inverse Fourier Transform
INVFT checks size and calls FFTI if size is a power of 2, else calls DFT.
input: level 2; vector
output:level 1; transformed vector
INVFT [D0C5]
<< FAST DUP SIZE 1 GET
IF LN 2 LN / FP
THEN -1 DFT
ELSE FFTI
END
ROUN
>>
--------------------------------------------------------------------------------
Periodic Convolution
NOTE: Input vectors MUST be of same size!
input: level 2; vector
input: level 2; vector
output:level 1; convolved vector of same size
PCNV [AA10]
<< FT SWAP FT -> a b
<< a SIZE 1 GET -> s
<< 1 s FOR n a n GET b n GET * NEXT s ->ARRY >>
>>
INVFT
>>
--------------------------------------------------------------------------------
Aperiodic Convolution
input: level 2; vector
input: level 1; vector
output:level 1; convolved vector of size size(v1) + size(v2) - 1
APCNV [A2D]
<< -> a b
<< a SIZE 1 GET
b SIZE 1 GET + 1 - 1 ->LIST -> s
<< a s RDM
b s RDM >>
>>
PCNV
>>
--------------------------------------------------------------------------------
Two Dimensional Fourier Transform
NOTE: Number of rows, and number of columns, must be an even power of two.
Pad matrix with zeroes if necessary.
input: level 2; matrix
output:level 1; transformed matrix
FFT2[9C46]
<< DUP SIZE LIST->
DROP -> a nr nc
<< "[" 1 nr
FOR r 1 nc
FOR c 'a(r,c)' ->NUM
NEXT
nc ->ARRY FT ->STR +
NEXT
STR-> 'a' STO "[" 1 nc
FOR c 1 nr
FOR r 'a(r,c)' ->NUM
NEXT
nr ->ARRY FT ->STR +
NEXT
STR-> TRN CONJ
>>
>>
--------------------------------------------------------------------------------
Inverse Two Dimensional Fourier Transform
NOTE: Number of rows, and number of columns, must be an even power of two.
Pad matrix with zeroes if necessary.
input: level 2; matrix
output:level 1; transformed matrix
IFFT2[F08F]
<< DUP SIZE LIST-> DROP
-> a nr nc
<< "[" 1 nc
FOR c 1 nr
FOR r 'a(r,c)' ->NUM
NEXT
nr ->ARRY INVFT ->STR +
NEXT
STR-> TRN CONJ 'a' STO "[" 1 nr
FOR r 1 nc
FOR c 'a(r,c)' ->NUM
NEXT
nc ->ARRY INVFT ->STR +
NEXT
STR->
>>
>>
--------------------------------------------------------------------------------
Axis Shift for Two Dimensional Fourier Transforms
Shifts axis from upper right corner to center of matrix, by swapping quadrant
2 with quadrant 4, and swapping quadrant 1 with quadrant 4.
NOTE: Number of rows, and number of columns, must be an even.
Pad matrix with zeroes if necessary.
input: level 2; matrix
output:level 1; shifted matrix
FTSH[5EFD]
<< DUP SIZE DUP LIST-> DROP
-> a s rt ct
<< rt 2 /
DUP 1 +
ct 2 /
DUP 1 +
-> r2 r3 c2 c1
<< r3 rt
FOR r c1 ct
FOR c 'a(r,c)' ->NUM
NEXT 1 c2
FOR c 'a(r,c)' ->NUM
NEXT
NEXT 1 r2
FOR r c1 ct
FOR c 'a(r,c)' ->NUM
NEXT 1 c2
FOR c 'a(r,c)' ->NUM
NEXT
NEXT s ->ARRY
>>
>>
>>
--
Robert Yoder 306 Hawkins Graduate House, West Lafayette, IN 47906 (317)495-6845
Internet: ryoder@ecn.purdue.edu "Flame all you want, We'll post more."
UUCP: ryoder!pur-ee Apologies to Jay Leno & Doritos.
Bitnet: ryoder@ecn.purdue.edu@purccvm.bitnet