Thread View: comp.dsp
1 messages
1 total messages
Started by Jacob Hirbawi
Fri, 09 Oct 1992 22:00
Re: FFT/DFT progs wanted...
Author: Jacob Hirbawi
Date: Fri, 09 Oct 1992 22:00
Date: Fri, 09 Oct 1992 22:00
129 lines
4629 bytes
4629 bytes
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 -------------------------------------------------------------------
Thread Navigation
This is a paginated view of messages in the thread with full content displayed inline.
Messages are displayed in chronological order, with the original post highlighted in green.
Use pagination controls to navigate through all messages in large threads.
Back to All Threads