Several FORTRAN program modules for performing one-dimensional and two-dimensional discrete Fourier transforms, Mellin, and Fourier-Bessel transforms are described along with programs that realize the algebra of high speed Fourier transforms on a computer. The programs can perform numerical harmonic analysis of functions, synthesize complex optical filters on a computer, and model holographic image processing methods

Debabov, A. S.

Ikhabisimov, D. K.

Kolosov, B. I.

Usikov, D. A.

English

NASA Technical Reports Server

5
 
NASA TM-75678
NASA TECHNICAL MEMORANDUM 

B80-18749
(IASA-TM-75678) PROGRAMS FOR HIGH-SPEED 

FOURIER, MELLIN AND FOURIER-BESSEL
 
TRANSFORMS (National Aeronautics and Space
 
Unclas
Administration) 19 p HC A02/HY A01 CSCL.09B 

_G3/61 47190
 
PROGRAMS FOR HIGH-SPEED FOURIER, MELLIN
 
AND FOURIER-BESSEL TRANSFORMS
 
DtKjoTkh4hjsi mqv, A. S. Debabov, B. I. Kolosov, 
D. A. Usikov
 
Translation of "Programmy Bystrykh Preobrazovaniy Fur'ye,
 
Mellina i Fur'ye-Besselya, Academy of Sciences USSR,
 
Institute of Space Research,Moscow, Report Pr-418,
 
1978, pp 1-17
 
NATIONAL AERONAUTICS AND SPACE ADMINISTRATION
 
Washington, D. C. 20546 September 1979
 
REPRODUCED BY 
TECHNICALNATIONAL 
INFORMATION SERVICE 
US. DEPARIMEIT OFCOMMERCE VA.22161SPRINGFIELD, 
https://ntrs.nasa.gov/search.jsp?R=19800010472 2020-03-21T20:04:46+00:00Z
NOTICE 
THIS DOCUM.ENT HAS BEEN REPRODUCED 
FROM THE BEST COPY FURNISHED US BY 
THE SPONSORING AGENCY. ALTHOUGH IT 
IS RECOGNIZED THAT CERTAIN PORTIONS 
ARE ILLEGIBLE, IT IS BEING RELEASED 
IN THE INTEREST OF MAKING AVAILABLE 
AS MUCH INFORMATION AS POSSIBLE. 
STANDARD TITLE PAGE 
1. 	 Report No. _2. Government Accession No. 3. Recipient's Catalog No.
 
NASA TI14-75678
 
4. Title and Subtitle 5. Report Date
 
PROGRAMS FOR HIGH-SPEED FOURIER; MELLIN AND September 1979
 
AND FOURIER-BESSEL TRANSFORMS 
 6. Performing Organiotio Code 
7. 	 Author(s) 8. Ps:forming Organization Report No. 
D. K. Tkhabisimov, A.- S. Debabov, B. I. 10. WorkUnitNo.
 
-- Kolosov-D A. -Usikov­
9. 	Performing Orgno zation Name ond Address 11. Contract or Grant No,
 
-SCITRAN NASW-3198
 
Box 5456 	 13. Type of Report and Period Covered 
Santa Barbara, CA 93108 	 Translation
 
12. 	 Sponsoring Agency Name and Address' 
National Aeronautics and Space Administration 14. Sponsoring Agency Coda
 
Washington, D.C. 20546
 
15. 	 Supplementary Notes Translation of 
tProgrammy Bystrykh P.reobrazovaniy Fur'ye, Mellina i Furtye-Besselya,.
 
Academy of Sciences USSR, Institute of Space Research, Moscow, Repoit
 
Pr-418 1978, pp. 1-17
 
16. 	 Abstract 
We present a descriptioh of program modules for performing one-dinensiona
 
and two-dimensional discrete Fourier transforms, Mellin and Fourier-

Bessel transforms, and also programd for realizing the algebra of high­
speed Fourier transforms on a computer. The programs developed can be
 
used to perform numbdri&al h&rmonic .analysis.of functions, t6 synthesize

complex optical filters on a computer, modelling the holographic methods
 
of processing images;
 
Th programs are written in FORTRAN and are included in the library of
 
''
progr'anis of 'the -ideor images processing unit'"SOF in the Institute of:
 
Space Sciences of the USSR Academy of Sciences.
 
17. 	 Key Words (Selocted by Author(s)) Ia. Distribulion Statement 
Unclassified - Unlimited 
19. 	 Securit Cbossil. (of this ireport 120. Security Clossif. (of this page) 21. No. of Pages 2 
Unclassified 	 Unclassified 
 17
 
iii
 
NAA.l 
TABLE OF CONTENTS
 
INTRODUCTION ....................... ........... 1
 
1. Program for the Higa-Speed Fourier Transform. ...... . . .. . . 3
 
1.1. The One-Dinensional Version .............. ........ 3
 
1.2. The Two-Dimensional Version ............... ..... 4
 
2. The Fourier-Bessel Transform . . , t ...... 4
 
3, Cliculating the Bessel and Gamma Functions .,.. .. . 5
 
4, Auxiliary Programs .......... ..... ......... 6
 
CONCLUSION ....................... ...............7
 
REFERENCES .... ........ 8
 
APPENDIX ........ ...... 9......
 
Preceding page blank
 
PROGRAMS FOR HIGH-SPEED FOURIER, MELLIN
 
AND FOURIER-BESSEL TRANSFORMS
 
D. K. Tkhabisimov, A. S. Debabov, B. I. Kolosov,
 
D. A. Usikov
 
We present a description of program modules for
 
performing one-dimensional and two-dimensional discrete
 
Fourier transforms, Mellin and Fourier-Bessel trans­
forms, and also programs for realizing the algebra of
 
high-speed Fourier transforms on a computer. The pro­
grams developed can be used to perform numerical harmonic
 
analysis of functions, to synthesize complex optical fil­
ters on a computer, modelling the holographic methods of
 
processing images.
 
The programs are written in FORTRAN and are in­
cluded in the library of programs of the video images 
processing unit "SOFI"' in the Institute of Space
 
Sciences of the USSR Acad---emy of "Sciences.
 
INTRODUCTION /3*
 
The need for numerical harmonic analysis arises when solving integral
 
equations of convolution type. Here all the functions are given in discrete
 
form. Let fK) be the file of the readings for the initial function 
which is nonzero on the interval of the real axis. A pair of discrete 
Fourier transforms (DFT) of the file J_-has the form:
 
K-0
 
and
 
fa6 Al Ct$ (2)
 
Numbers in the margin indicate pagination of original foreign text.
 
1 
where N is the dimension of the files f (n) and f (m). By Means of the trans­
forms (1) and (2) we can compute continuous Fourier transforms and the coeffi­
cients of the Fourier series of the initial function 1i]. The Mellin transform 
of the functionf(z)'is defined by the formula 12]: 
On making the substitution c=et' under the integral sign in (3), we obtain 
where ,-Z' is the operator of the continuous Fourier transform. From "(4) it 
is evident that the lHellin transform is carried out with the help of (1). The 
Fourier-Bessel function fC() has the form: 
,, f S4 z)fcx) / (5) /4 
a 
In Ill a procedure is shown for computing (5) with the help of the transforms 
Cil)-(3). Applying the Mellin transform C3) to (5), we obtain 
/(Pt >0 ) 
or 
C6)
 
where
 
IWTYZJC7) 
0The coefficient is computed at the points where r=e (~c'o ) 
The case n = 0 reduces to the case n = 1 if we note that 
where 
2 
To calculate (7) we can use the representation 123i
 
where ,L_ / When _17 t we make use of the recursion formula: 
Bessel's functions are computed by means of (1), since
 
In the present paper, programs are described which realize the transforms 
(1), (2), (3) and (5). A description is also given of programs for the algebra
 
of high-speed Fourier transforms, programs for calculating the Bessel and gamma
 
functions, etc. The text of the programs appears in the appendix.
 
1. Program for the High-Speed Fourier Transform
 
1.1 The One-Dimensional Version
 
The programs DQFT1 and RQFT1 are designed to carry out direct and inverse 
DFT of one-dimensional files described in the user's program as a complex. In 
accordance with the principles of organization of the 'complex "SOFI" J4], access 
to the programs is accomplished as follows: 
1) In the basic program the user defines
 
COM/IRZX t,
 
CO/oIAP Al T 111 f4
 
Illegible in original foreign text.
 
3 
where A is a working one-dimensional file with the DFT program to be processed 
2) The file to be processed is entered into the working file A and accesses
 
one of the programs of the high-speed Fourier transforms, e.g., CALL DQFTl;
 
3) The Fourier transform of the initial file is placed in-the file A; the
 
error in the calculation does not exceed 10- 8 when
 
1.2 The Two-Dimensional Version 
The direct and inverse DFT of two-dimensional files is performed with the
 
help of the programs DQFT2 and RQFT2. Here not every variable undergoes (1) and 
(2). In the basic program it is necessary to write: /6
 
COMI'1 A/ FICT'3yL6-) w/ ,,' 76(4,9,8 
where A, B are the one-dimensional working files with the DFT program to be 
processed (/4=2 ,=2 T..Moreover, analogous to Sec. 1.1, a description of 
the programs used in the high-speed Fourier transform programs can be found in 
[3].
 
When formulating the basic program, the user must designate the library
 
of the complex for which the card
 
must be inserted after the card
 
2. The Fourier-Bessel Transform
 
The Fourier-Bessel coefficients (5) are computed by the program 
F -FY4 L/g .47T )" Here F is the file for the readings of the 
initial function fCx) at the points X- e , NT is the number of readings on 
the interval 4 5-Z /IY_8 is the subscript of the Bessel function, W is tfre 
4 
file of the function , TM = tmax . For access to the program FTFB, 
the user must write in the basic program 
r¢tlP/VK F~m)ltt?rn) 
o-nroM. E _cT (3g_7, '," 
where m is the dimension of the file. The result is located in the file F.
 
The values of the function Wh(W) are calculated by means of the program
 
,Fr7-1z T(v TAzL& T . Here W is the desired file of values of the function 
_T is the number of values, A n..T ..... the function IMO) 
is computed at the points ITn, the program FLT, the 
gamma function is calculated with the help of series (cf. Sec. 3); however, it 
is also possible to make use of the representation (8). The programs FTFB and 
FILT are accessed by means of the operator CALL. 
3. Calculating the Bessel and Gamma Functions 
The program BESF(BF,TF,NF) is calculated by means of formula (9) NF of the 
Bessel functions at the point TF: 
..,W-T),47rF)
 
and the result is placed in the file BFCNF).
 
A.co.ute..h function 6 
The program the incomputes at
 
the points -%,AT -where M-7 and are the upper and
 
lower integration limits in (8). The result is located in the file G(NT). Be­
fore accessing programs BESF and GAlMA1, it is necessary to describe:
 
PEX )COI M qcy 
The program GAMA(Z) calculates the gamma function at the point ­
by means of the formula [5]:
 
/7 
5 
The result is located in Z.
 
4. Auxiliary Programs
 
Below we give descriptions of programs whose use simplifies access to the
 
high-speed Fourier transform programs and accelerates the process of formulating /8
 
a basic program. All programs work with complex files. 
1) hI)AE 4 X A14 A /9- the program selects line number K from file 
ACNX,NY) and when L = 0, copies the information from file A into file B(NY); when 
L = 1, it copies the information from file B(NY) into file A. 
AIX-- A2 the program selects column number K from 
file ACNX,NY); when L = 0, it copies the information from file A into file BCNX); 
when L i,1 it copies the information from file B(NIX) into file A. 
2) C//AIAjeA/x 
3) ou COM6') A1/V -- the program copies the ififormation from the 
COMON-block of element K in file C into file A(NXNY). 
4) TO CC/I £ -fQ/k1-K2 __ the program copies the information into the 
CONMON-block of element K in file C from file ANXNY). 
5) JfVZT4 6 P1W 'V) :- the program multiplies out the two two-dimen­
sional files A(NXNY and B(NXXNY), component by component, The result is located 
in th-e file S(NX,NY). 
-6 A'P( g8, ,Vr, A(S -- the program transposes the matrix A(NX,NY). 
The result is located in the file B(NX,NY).
 
Ff)>9OAXO(Y4-A " Y2-- the program forms the complex conjugate of 
the matrix A(NX,NY). The result is located in the file B(NX,NY),
 
178u.e0 V (., /// HY 'Ror R this is a program for the linear in­
terpolation of a function, given in Cartesian coordinates by the file of values
 
A(NX,NY) on a network with sampling steps HX and HY along the axes OX and OY,
 
respectively, into the polar exponential network (.e,y), where
 
' _r=" 4e aO Y- =J_C Oz9
 
r/e= e,t--- number of subdivisions of the network with respect /9CO 
6 
to the variable t, NF -- number of subdivisions of the network with respect to 
the variable f. The result is located in the file B(NT,NF). 
9) RYVO La>(F;,A/ V)XJ/Y/) is a program of linear interpolation 
of the values of a function given in the Cartesian file F(6di/) into the 
Cartesian network rotated through the angle AL relative to the origin of coordi­
nates. The origin of coordinates corresponds to the point 't'1 g ) 
NX and NY are the numbers of the nodes of the new network with respect to each
 
variable (the sampling step is not changed); V(NXNY) is the file of values of the 
initial function in the new network.
 
10) COR_ F( P* V_ _. . A/X, /V/s1. -- the program calculates the correlation 
function of the two images in the plane given by the files F(NxNY) and V(NX,NY) 
of the form: 
-~1X 1 jpc-)Wix) c/JX (10) 
' an.where 3 are vectors with components 2 d respectively. In 
the result, the file S(NX,NY) of the readings of the correlation function is ob­
tained.
 
CONCLUSION
 
In the present paper the programs described are included as object modules
 
in the complex "SOFI" and are intended for various problems in spectral processing;
 
recognition of images, refinement of images, the photographing of smog, the effects
 
of blurring, for various statistical and physical problems required when applying
 
the algorithms for the Fourier, Mellin and Fourier-Bessel transforms.
 
The authors express their deep appreciation to V. G. Zolotykhin for useful
 
discussions.
 
7
 
/10 REFERENCES 
1. D. A. Usikov. Prieneniye abstraktnogo gapmonicheskogo analiza dlya
 
bystrogo racpoznavaniya izobrazheniy (The Application of Abstract
 
Harmonic Analysis to High-Speed Identification of Images). Preprint
 
The Space Sciences Institute of the USSR Academy of Sciences, Pr-335,
 
Moscow, 1977.
 
2. 	N. Ya. Vilenkin. Spetsial'nyye funktsii i teoriya predstavleniy grupp
 
(Special Functions and the Theory of Groups).. Nauka, Moscow, 1965.
 
3. 	A. S. Debabov, D. A. Usikov. Spektral'nyye preobrazovaniya v komplekee 
obrabotki fotograficheskikh izobrazheniy "SOFI" iz YeS-EVM (Spectral 
Transforms in the Complex for Processing Images, "SOFI' on a YeS [Unified 
Series[ Computer), Part 1, Preprint, The Space Sciences Institute of the 
USSR Academy of Sciences, Pr-383, Moscow, 1977. 
4. 	V. G.-Zolotykhin, V. I. Kolosov, D. A. Usikov, et al. Opisaniye sistemy 
program natematicheskoy obrabotki na EVM Yedinoy Serii CYeS) fotografiches­
kikh izobrazheniy Zemli, poluchennykh s kosmicheskikh apparatov (kompleks 
"SOI") (Description of a System of Programs for Mathematically Processing 
on a Unified Series (YeS) Computer Photographic Images of the Earth Ob­
tained from Space Vehicles (The Complex "SOFI)). Preprint, The Space 
Sciences Institute of the USSR Academy of Sciences, Pr-336, Moscow 1977. 
5. 	Yanko, Ye., Emde F., AosfrF. Special'nyye funktsii (Special Functions).
 
Nauka Press, Moscow 1964.
 
8 
APPENDIX
 
DQFT 1 
Cl 
03 
C dHC 
kh-N 
IP!(394),NC,4,f,NlAlKCC( ) 
G405 
14 
C? 
Ca 
N2:412Ns35OZ-1 
CaVl fT' 
-%-4 
I ' ,th 
I? 
5 
16 
V9 
19 
14 
K:O 
CALL trio 
NC'i 
CALL FoTD 
NKC 
C(N(.1)(CCN 
96to C(SK. ,*)-C(NA4l)C NAh*ZIC(kxh s.%,)CC 4AX)C(HA*flhl 
:Z 
t. 'a 
CC.-
IETURN 
C.(1)- NA. 'I -C(NA ) 
21 C(hXA).C(MAt) Z 
104LEOUTI1 1 0FT1I 
13 
1& ~C(A:y0.(C.
3C(%A*I'I)aC.(CIv. 1)CC' 12-41)3.D)C(NK2.*1) 
C(-A.1).P *(C(%K2.*1I).CIUKaNI.I) 
s8 %cZrnC 
V CALL tI I 
13 CALL trio 
5 asc ? TI % 
a6 C( *K*.j-P).C(NA.t-1) 
ORIGINAL PAGE IS
 
JF POOR QUALITY
 
flOFT 2 
Cl~ 04o. FTCX( 304) ,fC.,flfl.,NA,fls,C(iI 
Z6 CALL I'kECO 
C9
10 
1A 
15 
17 
ii 
CALL T'.ECO 
4 .A 2 4 -
I.C.C. 
CALL uPTD 
4'.CALL FFI0 
X.eI ,ILOAT(hk? 
I C l ' ( ) x 
2 
14 
CLL 
CALL 
IhAECO 
P4ECO 
04 
O 
0010 
a 
C A L L 
CALL 
c -
CALL 
l s q c o 
?%AF(C 
Ci o1 
6 
17Is 
19 
4C.5 
CALL 14APCOSh-44s . 
~ 
DL P4 
O9 
S00Cs 
6 
r 
09 
10o. 
CrTurFX, 
1' .1 
9J)CL.EQ1j) 
60101 
2 A(C.)'t(J) 
I~u~ 
o 
FTFB 
01 5SLBROUIINE fTFK(F,llT,N8,v.T0.1XJ 
02 COPLF'--iP ~4 
0 5 ITTM-1O0 
06 OTsTT/FLOAT (?) 
0? NI.NI? 
Ca DO.IoI-,)'rk 
T)
()()EP
10 

11 C(2-*[ - tUEAL(fflt)l 
12 1C(2eZ)aAINAC(Pi)) 
I j CALL D0571 
14 OJ 2 t''Ia. 
is2 P(1)*C"DLY(t(2-3.j)CC2'1)3 
It 00 3 1 1,NI 
17 11 .NI 
I9a( F ( Itj s I I 
19 3 F(11);i(t1)MUI13
 
21C(alI-I)-REAL(FSCIU
 
24 CALL ft2FT 
RETURN4t7 
T 
F IL 
SUAROUTtINFjLCW4 NT,"B.TT)
 
COJVPLEX 
WCNIZ2
0! 

AOSALOO(2.)
04 
 p1k2.*3.i4l59265
o5 
 N. I
04 
 N1.fliz0? No-O'Ol(N2)
08 
 I4.(HB-NO)/Z-I
09 
 001T 'N 
0 1;(1-11?)TT
-IM 

14
 
is ilCMP LX tc5,041
 
I& CALL GA)MA(Z
 
I$ ~(I);
 
1 .EU.2) GOTO3
19 

pTj)'FLOAT(N3) /2.
 
22*Z 2.CPLX(XN.08)
 
2D1j 2 
4IW(l).Wfl)-2
 
-4 RR-C0S(2::OM:AO)
 
4tZtC 

2) 
RPI.X (RA flR I) 
ZIbCMPLX(FL0AT(ile),.2.0nM)
a? 

29 RR-REAL(W(T))
 
3D 1l'CMPLX(RRRl) 
3.
 
RETURN3$ 
ORIGINAL PAGE IS 
OF POOR QUALITYLI 
COMMO0N 1 1 1 3)I1 ( ) C l
 
06
 
bt Ttc1. , 
09 RRaCOS(1 F'?k Tx)) 
I.SC I-21N TY 
13 CAL 1 DOFTI 
14 00 2 I'~ h
 
17 FROI 
GAMNAI 
SUIRCUTZIE GAMMAI(G,NT,TO.T(,A8)
ClcOPt'A OtNI) A
 03COPMO~N C( ) nic.cl
 
D1 P .0)/IFILOAT (N)
 
a 0TC.1- I40
 
C? 

RS*T -XP(TX)
 
11 ZP(RAlL.165.) QOT02

12 OCC)to' 
Is CI)SCTM-To)-AP(R)
 
17 CALL. TOLOM(6,N,1,1)
 
it CALL oO;ri

18 CALL OUlCOH(q,N,1 .11
0 00 3 titNI
 
a1 A'.OC)
 
Z3 3 0 i)4 
1) 30RUTINE AMMA( A 
Y? rOPPLpx Z,s
C N-'KPo COMMO-(iPEAEA g flPObi3REAENII$ 
03 nw100o
 
04 RK'IALOII(FLOATCN)) 
05 R REAIIZ) 
RR CS R A) 
, CJ.SItJCfl(PA)
 
103o3'EXO. q. A) /2
 
11 'S4F Li(l)/CZ-FLOArtifl
 
Iu' 
EN 
01 su UIZ0ullECLHfCA t,NX,NYEKL
 
02 COPOLfX A(4XNY) .(x)
 
03 D0 1 lldAx
 
04 IF(L.tQ.l) 00102
 
05110"i I1 lx
cO
06 
.7 ? AK%,V)-b(j)
DO c h li
08 *COVTIqt8
 
10 fP4
 
su9tIoUl t AMUL1(ASNtNYJ
 
Li~COPtrXF A(N ,f),q(q7,HY).s(w'Nfl 
o0 tj1,Ihy034 DO I j TVA 
Cs S( .,)A(l.J)t(.l4J) 
06 R ltiI" 
07 END 
TACO"
 
ttP'OL T 14! T0t0"t,.NAVYtl
 
CC-PtF% Af(h,NY) 
J CCWUO% F~tT(119).'(t)
 
4 00 I jt1,6Y
 
as 00 1L 1
 
C L -fkL A(W. )) 
ou' C.P' 
C 9 
1 
c-i 
.3 
C ( L tI JA 1 AY 
l Ut' 
Lt'.. 
1 " | 
fie 
3 
e 
)' 
J ' ) 
r'Ii ,bCO 
s gLj'(t)0To(ANY 
.8 - (.til.C'tECCL CC. 
C fie 
f 1"R 4 j C 
-- AA if,,jiA,(At.Jf 
4 % 
ORIGINAL PAGE IS
 
OF POOR QUALITY
 
PoULR
 
01 $4ROUTINE PDULRCAIIXHY.NXNYBNTNFRoR)
03 COMPLEX A(NXpjY) ,( NTfFr 
XOO.(003 NX- CIX
0& YO.Stp-I)* j ATy

05 DT-ALOOa, RPI FLOAT(NT)
06 00 1 Lj t.NF
0? cO 
i rc ,N1 
S06 BCKL)O.
 
9 00 2 L1.h F
10 - O ri,IJX 
.11 RKo.OEXP(,K1).DT) 
F2P2.'3.14159265*(L-1)/FLOAT(UF)
4 Y RK C6 1 3
 
6 J +,N7 (Y-YO)/Hy)
 
IS ji-J*1
 
I CI..0.NX)G To2 
I 1.L GY) a
 
2 CONTIHU!
 
4- BOURN 
REVOL
 
01 8UBROnJTflE REVQLCF NO,N ,HYAI.I
 
1 PRO co(E
01CmhiFH'i,J RX 1+, .(jJ.~p
NYJ~j

suBJou f
I 
J0 .NY 1i )
1 =- lCAL4 O C 

FRO L9 
'" j OUT IRE CORF
.SUR CORF(F.Vl5,NXNYI
 
COMMON FIC ( 9.) .R,"N BZt Z~ C
 
CA L 70 OF,Hl~y l)
 
log CALL Va 'CI Y'! EN Y,'I)0Q CAULL
09 ALL AQ172 NNy I
 
i ! CALL Y.I ' V, 
13I a Ml C i NY,iJN)
 
I&1 FALt T
 
13A LObI¢C SNY
 
OFPo 40 
O'OOR , I 
V41~l~ 


Programs for high-speed Fourier, Mellin and Fourier-Bessel transforms

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server