From: jcbhrb@nic.cerf.net (Jacob Hirbawi)
Newsgroups: comp.dsp,comp.sys.transputer
Subject: Re: FFT/DFT progs wanted...
Date: 9 Oct 92 22:00:42 GMT
Followup-To: comp.dsp
Organization: CERFnet

In comp.dsp <1992Oct9.154544.20151@csd.uwe.ac.uk>
PJ Pycroft  <pj_pycro@csd.uwe.ac.uk> writes

> I'm in the progress of setting up a project on noise recognition systems,  
> recognize some form of noise such as siren (police, fire etc.), alarm
> signal or any other warning noise which will alert a deaf person.
> 
> I'm will be using Signal Processing and Neural Network for detecting
> above noise and will be using parallel processing such as Transputers.
> I won't be using DSP chip's cos' of the limitation of our facility
> but Transputers will do the job....

> I would be grateful if you could help me ... where can I find the FFT
> as well as DFT programs (128, 256, 512 points and perhaps more!!.)
> in 'C' source code and if there is any possibility for me FTP' these
> programs.

> And, any other info' on above subject would be greatly appreciate...

> I'm a profoundly deaf postgraduate student doing M.Sc. in Parallel Computing
> Systems.

> Thanks for your help.

There is a large number of C fft routines available over the net. The best
way to look for them is of course through archie which will probably give you
more sources than you care to find. Anyway, since you mention that you will be
using Transputers how about coding the fft in Occam! -- not that I prefer
Occam over C but for relatively compact signal processing code like fft's
there really is little difference between the two; now for more complex
software development it's another story ;-) . Anyway I ran the code listed
below on a T800 and timed it at about 80 msec for a 1024 complex array -- no
where near a legitimate DSP (< 4 msec) but still not too bad. It consists of
three parts: (1) generate the 'twiddle factors' in routine Omega;(2) do a
fast in-place bit reversal in routine BitReverse; and (3) do the fft operations.

To use it you need to declare four arrays: the real and imaginary part of your
data of size N; and the real and imaginary parts of the 'twiddle factors' of
size M where 2^M = N; for examples:

  [512]REAL32 DataReal,DataImaginary:
  [9]REAL32 OmegaReal,OmegaImaginary:
  SEQ

    DataReal      :=   ( define or acquire your real data )
    DataImaginary :=   ( define or acquire your imaginary data )

    Fft(DataReal,DataImaginary,OmegaReal,OmegaImaginary)

    -- DataReal and DataImaginary now correspond to the spectrum of the data.

Good Luck!

Jacob Hirbawi
JcbHrb@CERF.net

-- CUT HERE -------------------------------------------------------------------
PROC Omega([]REAL32 w1,w2)
  INT i1,j1 :
  REAL32 a1 :
  SEQ
    i1:=1
    SEQ j1 = 0 FOR SIZE w1
      SEQ
        a1:=3.1415926345 (REAL32) /(REAL32 TRUNC i1)
        w1[j1]:=  COS(a1)
        w2[j1]:= -SIN(a1)
        i1:=i1+i1
:
PROC BitReverse([]REAL32 x1,x2)
  INT i1,i2,i3,j1,j2,n :
  REAL32 t1,t2 :
  SEQ
    n := SIZE x1
    i2:=1
    SEQ i1 = 1 FOR (n-1)
      SEQ
        IF
          i1<i2
            SEQ
              t1:=x1[i1-1]
              t2:=x2[i1-1]
              x1[i1-1]:=x1[i2-1]
              x2[i1-1]:=x2[i2-1]
              x1[i2-1]:=t1
              x2[i2-1]:=t2
          TRUE
            SKIP
        i3 := n/2
        WHILE i3<i2
          SEQ
            i2:=i2-i3
            i3:=i3/2
        i2:=i2+i3
:

PROC Fft([]REAL32 x1,x2,w1,w2)
  INT i1,i2,i3,j1,j2,n,m,k1,k2:
  REAL32 u1,u2,q1,q2,t1,t2:
  SEQ
    Omega(w1,w2)
    BitReverse(x1,x2)
    m := SIZE w1
    n := SIZE x1
    i1:=1
    SEQ j1 = 0 FOR m
      SEQ
        i2:=0
        WHILE i2<n
          SEQ
            u1:= 1.0 (REAL32)
            u2:= 0.0 (REAL32)
            SEQ i3 = i2 FOR i1
              SEQ
                k2:=i3+i1
                q1:=(u1*x1[k2])-(u2*x2[k2])
                q2:=(u2*x1[k2])+(u1*x2[k2])
                x1[k2]:=x1[i3]-q1
                x2[k2]:=x2[i3]-q2
                x1[i3]:=x1[i3]+q1
                x2[i3]:=x2[i3]+q2
                t1:=u1
                t2:=u2
                u1:=(t1*w1[j1])-(t2*w2[j1])
                u2:=(t2*w1[j1])+(t1*w2[j1])
            i2:=i2+(i1+i1)
        i1:=i1+i1
:
-- CUT HERE -------------------------------------------------------------------