Conference rusure::math

The best reference I know of off hand is the IEEE book 'Programs for Digital
Signal Processing'... I don't know who has the Magtape with the programs
here, but they are definitely available.  One of the programs, a very
fast radix 4 routine was done by Robert Morris, here in corporate R+D,
but he no longer works for DEC.

- Jim

There appers to be an "old reliable" FFT in the Laboratory Subroutine Package
(LSP) put out by the LDP group. It will run on a MINC so should be runnable
on PRO-350's etc. I think it's a MACRO version. I now have a copy of FOUR1,
a small FFT program written in FORTRAN. I am told there is a FOUR2 and a
FOURT about somewhere, but I don't have them.

Lynn Yarbrough

There is a good article on FFT's, entitled 'FFT Algorithms for Vector
Computers', by P. N. Swarztrauber, including FORTRAN algorithms, in
"Parallel Computing" Journal #1 (1984), recently published by
North-Holland, and available in the ZK library. The Cooley-Tookey, Pease,
and Stockham Algorithms are given, with analysis of speed and storage
requirements on parallel/vector processors as well as serial processors. 

Lynn Yarbrough

    VAXLAB has a signal processing package called LSP which has real
    and complex forward and inverse FFT's. A copy of the software can
    be obtained from LDP::BELANGER. The routines are based on the IEEE
    algorithms with some VAX specific performance enhancements. Performance
    data can be obtained from LDP::SCHURE.
    
    Amr
    

From the June 1987 IEEE Computer Magazine (P. 82):

"	ALGORITHM REPLACES FOURIER TRANSFROM

" Ronald Bracewell, a professor of electrical engineering at Stanford 
University, has patented an algorithm to replace the Fourier transform.

" He named the algorithm after Ralph Hartley, because of Hartley's work at 
Bell Laboratories in 1942, but the Patent Office calls it the Bracewell 
Algorithm. According to Bracewell, the former name seems to have stuck
among engineers. 

" Bracewell assigned the patent to the Trustees of Leland Stanford Junior 
University, which will license out the algorithm. According to the 
university, once the overhead is paid for adminstering the patent, 
Bracewell will earn one-third of the income from the licenses.

" The inventor is reputedly trying to convince other engineers at Stanford 
to build a chip containing the algorithm, which would give the university a 
patent over both the algorithm and the chip."

    I believe that Bracewell has published a textbook on the Hartly
    transform a year or two ago, and also he has published at least one
    paper on it - check the IEEE index of authors for his name in the past
    year or two (I don't have a copy handy so can't comment on the details
    of the Hartly transform).

    I am surprised that it is being patented though.

    - Jim

    Long ago, LDP did a nifty FFT routine for the PDP-11 as part of
    a package called LSP-11 (for Laboratory Subroutine Package). The
    MACRO-11 source was translated to MACRO-32 by SWS/Japan sevral years
    ago. In the course of my own search for a VAX FFT routine, I dug
    up the following information on the status of that code:
    
    (Note: I had asked specifically about licensing issues, 
    since I was looking to get some code for a friend outside
    the company...)

From:	ROBUST::ROBUST::MRGATE::"RUGGED::AKOV00::TKOV00::TKOV50::DECMAIL::43675"  9-FEB-1988 12:28
To:	MRGATE::"LDP::FORBES"
Subj:	Re: SWS/Japan FFT routine conversion?

From:	NAME: MATSUNAGA
	INITLS: MASAMITSU
	FUNC: SWS PSS BUSINESS
	ADDR: TKO
	TEL: 03-989-7046 <43675@DECMAIL@TKOV50@TKO>

To:	NAME: FORBES
	INITLS:   <FORBES @LDP@MRGATE@ROBUST@RUGGED@AKOV00@TKOV00 @TKOV50>
	RONALDO FORESTI @AKO @TKOV50
	NAME: ITOH
	INITLS: NAOKI <43756 @DECMAIL @TKOV50 @*>
	NAME: SEKIGUCHI
	INITLS: HISAYOSHI <84990 @DECMAIL @TKOV50 @*>

	Regarding to FFT routine, we have converted LSP-11 (not 
REAL-11) to make it run on VAXes. LSP-11 contains FFT routine.

If you would like to get this program, we can offer this without 
charge provided that:

1. You use the original version of LSP-11 document in English. Because 
the functions are not changed, you can use this without problems.
(I.E. We can not provide the documentation.)
2. You manage licensing issue by yourself. You don't ought to pay the 
license to us. However, that is Digitals property. So, please protect 
our data right.
3. You deliver the support service, if necessary. We can not deliver 
any support service for this product. We believe that this product is 
very stable and you can look into the source program to analize 
problems even if it happened.


	If you can accept the above and need the program, please 
inform me. Regarding to delivery method, the network access will be 
easiest way.

	Best Regards,
	Matsunaga Masamitsu

    Re .1:
    
>here, but they are definitely available.  One of the programs, a very
>fast radix 4 routine was done by Robert Morris, here in corporate R+D,
>but he no longer works for DEC.  ^^^^^^^^^^^^^
				  ^^^^^^^^^^^^^
    
    Hah!  I have been wondering when and where he worked for us!

    I obtained a copy of those routines some time ago, so if someone has
    the need let me know and I can place them on line.

    I also have Marple's recent took on spectrum analysis, which has a
    useful diskette of FORTRAN programs.

    Finally, you should know that the great speed gains are in going from
    the O(N^2) algorithm to *any* FFT routine - the fancy knotted code
    routines are only a certain factor faster; most people could just
    use a generic FFT from Numerical Recipes, or a DSP textbook like
    Oppenheim and Shaeffer, and get useful results quite quickly.  So if
    there is trouble to obtain the fancy FFT's - try a simple minded one
    first!

    - Jim

The stuff the DXML group has been doing recently is quite impressive. For 
example, from the DXML conference:

	"A second type of FFT is in the process of being released as
	part of DXML which is known as the group FFT.  This can be used
	when the user wishes to perform the same length FFTs across
	several data sets at the same time, (ie A(N,1024) can use the
	group FFT to compute all N FFTs of length 1024.  The advantage
	shows up in the performance.  Here is some Preliminary Data on
	the VAXvector 9000 using the group FFT;

	Complex-Complex D.P. 1024 point in 1.29 msec gives 39.8 MFLOPs.

	The person to contact in DXML is Roy Ho who has done a great
	job in improving the FFT performance using this new group
	technique.  NOTE: These are all double precision results,
	if single precision is desired, the performance does improve
	on the VAXvector 9000.  For the VAXvector 6000 machine single
	precision is almost twice as fast as double precision.  I hope
	this data helps."

    re .10,
    
    Well, as of last Monday, I have joined the DXML group.  I have
    forwarded a copy of .0 to him.  Glad to hear a long time friend (we
    used to work in the High Performance Technology group) receive this 
    wide recognition.
    
    Eugene

    I'm having a problem with the FFT algorithm from NUMERICAL RECIPES
    (fortran,pascal,c, they're all the same). It seems that the imaginary
    part of the result is the wrong sign.
    
    I'm generating a periodic time domain waveform, transforming to freq
    domain, applying a filter function to the components, and transforming
    back to time domain. The problem is that the result appears to have a
    negative time shift instead of the expected positive shift. So, I took
    my time domain source waveform and used Excel to calculate the FFT.
    Comparing the imaginary part of Excel's transform and Numerical
    Recipes', the sign of each component is opposite. 
    
    Has anyone else noticed this? I can't believe that an algorithm this
    old could be wrong and not have been caught by now. (I first used the 
    algorithm from the first edition of "Numerical Recipes" and translated it
    myself into C, then looked it up in the 2nd Edition of "N.R. in C"
    identical routine) Could I be interpreting the result wrong somehow?
    
    Thanks,
    
    Steve M.

    Can't help you directly, but it seems to me that I recently saw a
    title something like, "Error in Numerical Recipies FFT?" in one of
    the USENET math groups (probably, sci.math or sci.math.numeric).  I'll
    check tomorrow to see if it was recently enough to still be active.

    Do you have the supplement?  This is a volume of test cases and
    sample calling sequences.  I haven't seen it for the second edition,
    but I have one (Pascal) for the first.  I can see what they have
    for FFT if you like.

					Topher

I had the same problem with the FFT in Numerical Recipes.  Basically, it's because
their background is in physics, where the assumption is that a frequency domain value has
an implied e^(-i\omega t), rather than the electrical engineering e^(j\omega t).  Just
flip the sign of the imaginary part on the transform (or on your filter!) and you'll
be all set.
					-Michael

    I had the same problem with the FFT in Numerical Recipes.  Basically,
    it's because their background is in physics, where the assumption is
    that a frequency domain value has an implied e^(-i\omega t), rather
    than the electrical engineering e^(j\omega t).  Just flip the sign of
    the imaginary part on the transform (or on your filter!) and you'll be
    all set.
    
					-Michael