home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #27 / NN_1992_27.iso / spool / comp / sys / transput / 1214 < prev    next >
Encoding:
Internet Message Format  |  1992-11-18  |  3.1 KB
  1. Path: sparky!uunet!ornl!rsg1.er.usgs.gov!darwin.sura.net!zaphod.mps.ohio-state.edu!sdd.hp.com!usc!news!nic.cerf.net!jcbhrb
  2. From: jcbhrb@nic.cerf.net (Jacob Hirbawi)
  3. Newsgroups: comp.sys.transputer
  4. Subject: RE: FFT OCCAM CODE WANTED !!!
  5. Message-ID: <3779@news.cerf.net>
  6. Date: 19 Nov 92 03:01:52 GMT
  7. Sender: news@news.cerf.net
  8. Organization: CERFnet
  9. Lines: 105
  10. Nntp-Posting-Host: nic.cerf.net
  11.  
  12. In comp.sys.transputer <1992Nov18.173159.5761@rmece49.upr.clu.edu>
  13. Roque Pagan <rpagan@rmece02.upr.clu.edu> writes:
  14.  
  15. >    I'm trying to locate some FFT source code in OCCAM.  Any
  16. >    information would be greatly appreciated.
  17.  
  18. This is from a response I posted a few weeks ago:
  19.  
  20. ... The code consists of three parts: (1) generate the 'twiddle factors' in 
  21. routine Omega,(2) do a fast in-place bit reversal in routine BitReverse,
  22. and (3) do the fft operations.
  23.  
  24. To use it you need to declare four arrays: the real and imaginary part of your
  25. data of size N; and the real and imaginary parts of the 'twiddle factors' of
  26. size M where 2^M = N; for examples:
  27.  
  28.   [512]REAL32 DataReal,DataImaginary:
  29.   [9]REAL32 OmegaReal,OmegaImaginary:
  30.   SEQ
  31.  
  32.     DataReal      :=   ( define or acquire your real data )
  33.     DataImaginary :=   ( define or acquire your imaginary data )
  34.  
  35.     Fft(DataReal,DataImaginary,OmegaReal,OmegaImaginary)
  36.  
  37.     -- DataReal and DataImaginary now correspond to the spectrum of the data.
  38.  
  39. Good Luck!
  40.  
  41. Jacob Hirbawi
  42. JcbHrb@CERF.net
  43.  
  44. -- CUT HERE -------------------------------------------------------------------
  45. PROC Omega([]REAL32 w1,w2)
  46.   INT i1,j1 :
  47.   REAL32 a1 :
  48.   SEQ
  49.     i1:=1
  50.     SEQ j1 = 0 FOR SIZE w1
  51.       SEQ
  52.         a1:=3.1415926345 (REAL32) /(REAL32 TRUNC i1)
  53.         w1[j1]:=  COS(a1)
  54.         w2[j1]:= -SIN(a1)
  55.         i1:=i1+i1
  56. :
  57. PROC BitReverse([]REAL32 x1,x2)
  58.   INT i1,i2,i3,j1,j2,n :
  59.   REAL32 t1,t2 :
  60.   SEQ
  61.     n := SIZE x1
  62.     i2:=1
  63.     SEQ i1 = 1 FOR (n-1)
  64.       SEQ
  65.         IF
  66.           i1<i2
  67.             SEQ
  68.               t1:=x1[i1-1]
  69.               t2:=x2[i1-1]
  70.               x1[i1-1]:=x1[i2-1]
  71.               x2[i1-1]:=x2[i2-1]
  72.               x1[i2-1]:=t1
  73.               x2[i2-1]:=t2
  74.           TRUE
  75.             SKIP
  76.         i3 := n/2
  77.         WHILE i3<i2
  78.           SEQ
  79.             i2:=i2-i3
  80.             i3:=i3/2
  81.         i2:=i2+i3
  82. :
  83.  
  84. PROC Fft([]REAL32 x1,x2,w1,w2)
  85.   INT i1,i2,i3,j1,j2,n,m,k1,k2:
  86.   REAL32 u1,u2,q1,q2,t1,t2:
  87.   SEQ
  88.     Omega(w1,w2)
  89.     BitReverse(x1,x2)
  90.     m := SIZE w1
  91.     n := SIZE x1
  92.     i1:=1
  93.     SEQ j1 = 0 FOR m
  94.       SEQ
  95.         i2:=0
  96.         WHILE i2<n
  97.           SEQ
  98.             u1:= 1.0 (REAL32)
  99.             u2:= 0.0 (REAL32)
  100.             SEQ i3 = i2 FOR i1
  101.               SEQ
  102.                 k2:=i3+i1
  103.                 q1:=(u1*x1[k2])-(u2*x2[k2])
  104.                 q2:=(u2*x1[k2])+(u1*x2[k2])
  105.                 x1[k2]:=x1[i3]-q1
  106.                 x2[k2]:=x2[i3]-q2
  107.                 x1[i3]:=x1[i3]+q1
  108.                 x2[i3]:=x2[i3]+q2
  109.                 t1:=u1
  110.                 t2:=u2
  111.                 u1:=(t1*w1[j1])-(t2*w2[j1])
  112.                 u2:=(t2*w1[j1])+(t1*w2[j1])
  113.             i2:=i2+(i1+i1)
  114.         i1:=i1+i1
  115. :
  116. -- CUT HERE -------------------------------------------------------------------
  117.