Parallax Forums

Parallax Forums > Public Forums > Propeller Chip > Fast Hartley Transform

Forum Quick Jump

[ << Previous Thread | Next Thread >> ] | Show Newest Post First ]

ngl
Registered Member

Date Joined Sep 2006
Total Posts : 17

Posted 10/29/2009 5:22 PM (GMT -8)

There has been some discussion of the Fast Fourier Transform implementation on the propeller in another thread on this forum. One approach was detailed in the Propeller Wiki. I would suggest that a more practical approach to fft implementation is the Fast Hartley Transform (FHT), which works with real numbers, therefore running faster than the FFT. Furthermore, the FHT is its own inverse function and the FFT can be calculated from the FHT. Ullmann (sepwww.stanford.edu/theses/sep38/38_29_abs ) has described a FHT algorithm (expressed in Ratfor) which is easily translated into Spin.

The attached spin file is a first implementation of the 1-dimensional HFT, which also includes methods for calculating the Power Spectrum and the FFT. ViewPort was used to debug the program and view the results. Things to do include providing a graphical data display, converting the subroutines (eg, butterfly) to pasm, timing studies with larger data sets, and implementing the 2-dimensional FHT.

Comments and suggestions are welcome.

Nick

Corrected spin file and url.

Post Edited (ngl) : 10/31/2009 4:49:28 AM GMT

File Attachment :
FHT.spin 5KB (application/octet-stream)

This file has been downloaded 17 time(s).

Kye
Kwabena Agyeman

Date Joined Feb 2008
Total Posts : 619

Posted 10/29/2009 8:05 PM (GMT -8)

Thank you for contributing.

Nyamekye,

Drone
Registered Member

Date Joined Mar 2007
Total Posts : 359

Posted 10/30/2009 4:04 AM (GMT -8)

Thanks for the work Nick. The link in the Ullman reference in your post and in the source spawns a 404 Not Found. How useful did you find the Bracewell book on FHT? I've seen this book for sale on the Web but wondered if it was too dated. I think Bracewell developed the DHT, but I'm wondering if the book is too old to cover FHT variants.

Regards, David

Post Edited (Drone) : 10/30/2009 12:13:52 PM GMT

hover1
I have 3 Propellers

Date Joined Sep 2007
Total Posts : 371

Posted 10/30/2009 4:14 AM (GMT -8)

http://sepwww.stanford.edu/theses/sep38/38_29_abs.html

Jamesx
Registered Member

Date Joined Aug 2005
Total Posts : 123

Posted 10/30/2009 5:17 AM (GMT -8)

Thanks, Nick.

I am very interested in this work. My problem is that I don't know much about FFT or FHT, but having this code and references gives me somewhere to start. I am looking forward to your successes at getting it into ASM, as well as tests of timing.

Jim C.

Hanno
Registered Member

Date Joined May 2007
Total Posts : 447

Posted 10/30/2009 11:41 AM (GMT -8)

Great job Nick!
(And nice use of ViewPort!) Got any screenshots? Any chance you could apply this to my "Speech Recognition" challenge?
Hanno

Co-author of the official Propeller Guide- available at Amazon
Developer of ViewPort, the premier visual debugger for the Propeller (read the review here, thread here),
12Blocks, the block-based programming environment (thread here)
and PropScope, the multi-function USB oscilloscope/function generator/logic analyzer

ngl
Registered Member

Date Joined Sep 2006
Total Posts : 17

Posted 10/30/2009 8:43 PM (GMT -8)

David: I find Bracewell very useful - I am still using it.

Hanno: ViewPort allowed me to quickly find programming errors. A great product. I will have to look at the speech recognition challange.
I would like to use the xy widget in ViewPort to display the FHT results. Any hints on how to do this?

There is a sign error in the FFT method. The last statement should be Im[k] := (-fx[k] + fx[len-k+2])/2

I have converted the program to a true object which also includes convolution & correlation methods. I will post as soon as I have a working demo.

Hanno
Registered Member

Date Joined May 2007
Total Posts : 447

Posted 10/30/2009 11:08 PM (GMT -8)

Hi Nick,
Thanks! The xy widget takes the values of two traces at the current timescale and plots them in xy fashion. Not sure what you want to do, but typically you would share at least 2 variables, select the proper timescale, and plot 2 channels.
Hanno

BR
Omelette du fromage!

Date Joined Dec 2008
Total Posts : 6

Posted 10/31/2009 8:12 AM (GMT -8)

Nick,

This is very cool. I don't have a copy of viewport (yet), so I made a version that uses the Parallax Serial Terminal. I set up the output so it can be directly imported into PLX_DAQ for easy plot/visualization. I also added sin, cos, step, and impulse functions to the examples you had so I could see the output with various waveform inputs.

I included a copy of my PLX_DAQ excel file in the archive...for anybody who'd rather get a fresh copy of PLX_DAQ, it is available here:
http://www.parallax.com/Portals/0/Downloads/sw/plx_daq_install.zip

The output from the various input waveforms isn't exactly what I was expecting...but that may be as much due to my ignorance of the subject area as anything else. Between you, Jay, and lonesock, and others, it's been a veritable frequency transform cornucopia here in the forums lately.

Looking forward to your next update.

BR

EDIT: updated archive attachment:
-Added a random (white) noise function, pulse function
-Added a simple measurement of #clock cycles required to complete H-xform
-Added magnitude and phase plots in excel PLX_DAQ worksheet
-Misc. tweaks

Post Edited (BR) : 11/2/2009 3:19:54 AM GMT

File Attachment :
FHT - Archive [Date 2009.11.01 Time 22.15].zip 84KB (application/x-zip-compressed)

This file has been downloaded 5 time(s).

ngl
Registered Member

Date Joined Sep 2006
Total Posts : 17

Posted 11/2/2009 5:22 PM (GMT -8)

Attached find a 2nd version of the FHT project reformatted as a basic FHT object and a demo program to exercise the FHT object.

BR: Thanks for your interest and additions to FHT. I have incorporated them, with modifications to accomodate the new FHT object, in the attached demo program. That is, your earlier additions. I look forward to making further additions to the demo program.

I have replaced attached programs with new versions containg BR's latest additions.

Nick

Post Edited (ngl) : 11/3/2009 2:07:22 AM GMT

File Attachment :
FHT_20091101.spin 4KB (application/octet-stream)

This file has been downloaded 15 time(s).

File Attachment :
FHT_Demo_1.spin 15KB (application/octet-stream)

This file has been downloaded 11 time(s).

ericball
Registered Member

Date Joined May 2007
Total Posts : 277

Posted 11/2/2009 6:36 PM (GMT -8)

I'm reading through US Patent 4,646,256 for the (Fast) Discrete Bracewell (Hartley) Transform and I'm not sure the RATFOR code is correct, particularly the initial bit reversal step. It might be a good idea to test it out by calculating the DHT (i.e. use the full O(N^2) calculation.), just watch out for scaling differences.

None the less, it looks like an interesting alternative to the typical FFT.

Composite NTSC sprite driver: Forum

NTSC & PAL driver templates: ObEx Forum

OnePinTVText driver: ObEx Forum

ngl
Registered Member

Date Joined Sep 2006
Total Posts : 17

Posted 11/2/2009 9:38 PM (GMT -8)

ericball: I did make one change to the bit reversal subroutine - butdis was initialized to 1 rather than 0 in Ullmann's document.
The code works in that I could reproduce all Bracewell's examples that I tried.

Nick

I also cross-checked the FFT methgod using the Descrete Fourier Transform in Mathematica.

Post Edited (ngl) : 11/3/2009 6:22:52 AM GMT

ericball
Registered Member

Date Joined May 2007
Total Posts : 277

Posted 11/4/2009 6:13 AM (GMT -8)

I've gone through Bracewell's patent and I think Ullman's algorithm must be a second generation algorithm. The Bracewell algorithm uses a 2 element butterfly, i.e. a1 = a0 + b0 * cas b1 = a0 - b0 * cas, rather than the 3 (4) element butterfly used by Ullman. I'm contemplating implementing the Bracewell algorithm in SPIN, but I'd like to know how many samples & sample width I should be targetting.

I'm hoping I can squeeze the algorithm into 240 instructions, leaving 256 registers for the element array, which would speed things up slightly (and make the bit reversal permutation trivial as it would happen on the load from HUB to COG). Moving the data out to HUB RAM means the array could be any size, but the elements would have to be LONGs. (Probably 16+16 signed fixed point.)

And what sample data do people use? 8 bit signed or unsigned, 16 bit signed, or some fixed point? The output is always signed.

Composite NTSC sprite driver: Forum

NTSC & PAL driver templates: ObEx Forum

OnePinTVText driver: ObEx Forum

BR
Omelette du fromage!

Date Joined Dec 2008
Total Posts : 6

Posted 11/4/2009 8:27 PM (GMT -8)

@Nick: thanks for posting your updated object/demo. I've been trying it out, seems to work fine. I took the liberty of making another tweak...I replaced all the "/2" and "/scale" operations with shift arithmetic right operations, also changed the scale to 1024 to facilitate this. Was worth about 8% improvement in execution speed, at the cost of code readability. Maybe worth the trouble, if for no other reason than for making the spin code match what would likely be the preferred implementation if you decide to covert it to ASM. It would be interesting to see how an ASM implementation of the FHT stacks up against Ale's FFT on the wiki page. The wiki claims 1024 sample FFTs @ 16 fps with a single cog.

@ericball: as to sample data, 16 bit signed seems to be a good fit for the prop. cessnapilot has been proponing a "Qs15_16" fixed point format (http://forums.parallax.com/forums/default.aspx?f=25&m=374493). I've been pondering what it would take to implement this for a digital filter project I've been experimenting with. Might be worth considering for your Bracewell algorithm, assuming you decide to put the data in the hub ram.

File Attachment :
FHT_20091101_sar.spin 9KB (application/octet-stream)

This file has been downloaded 5 time(s).

ericball
Registered Member

Date Joined May 2007
Total Posts : 277

Posted 11/5/2009 6:07 AM (GMT -8)

Yeah, I've come to the conclusion that storing the table in HUB RAM is the way to go, and not just so the table can be larger. I was thinking that storing the table in COG RAM would allow the HUB RAM table to be words or even bytes, but then there is a significant loss of precision when writing the output values back to HUB RAM; not good if you want to reverse the transform. So it's best to use the same format (probably Qs15_16) for both intermediate and final values.

Using HUB RAM also opens the possibility of having multiple cogs process the data simultaneously. There would need to be some mutual synchronization to ensure each stage is completed before the cogs start the next stage. But other than that, the algorithm parallelizes quite nicely.

Note: the scaling should be 1/sqrt(N) so the transform is it's own inverse. In my algorithm I'm going to shift the result as a part of every other stage, which should help with overflow.

Composite NTSC sprite driver: Forum

NTSC & PAL driver templates: ObEx Forum

OnePinTVText driver: ObEx Forum

ericball
Registered Member

Date Joined May 2007
Total Posts : 277

Posted 11/6/2009 5:43 PM (GMT -8)

BR said...
It would be interesting to see how an ASM implementation of the FHT stacks up against Ale's FFT on the wiki page. The wiki claims 1024 sample FFTs @ 16 fps with a single cog.

I've done a quick cycle count of my untested PASM FHT - 266 cycles * N * log₂N , or about 30 1024 sample FHTs per second at 80 MHz using a single cog. Now, the output is a Hartley Transform cas frequency domain, so if you want the Fourier eⁱ frequency domain you'll have to convert. Or you can just do all of the work in the cas domain directly. All the usual stuff should be possible, assuming you're dealing with real signals in the time domain.

edit 2009-11-12
It turns out the FHT starts doing the 4 element butterfly after two rounds of 2 element add/sub transforms. I didn't realize it until I started executing the third round by hand and discovering it didn't match the expected results. Then I went back through the patent and paid more attention to the initial theory section. Unfortunately, this makes the FHT more complex, the inner loop is now 364 cycles per sample, so around 21.5 1024 sample FHTs per second. (That's a minimum rate, actual may be slightly faster.)

I also discovered a bunch of problems with my math routines - signed values are a little trickier to handle than unsigned values.

I'm editing this to note my errors & progress (hopefully) without bumping the topic to the top.

Composite NTSC sprite driver: Forum

NTSC & PAL driver templates: ObEx Forum

OnePinTVText driver: ObEx Forum

Post Edited (ericball) : 11/13/2009 1:10:34 AM GMT

[ << Previous Thread | Next Thread >> ]

Forum Information

Currently it is Friday, November 20, 2009 11:06 PM (GMT -8)
There are a total of 393,737 posts in 55,521 threads.
In the last 3 days there were 82 new threads and 701 reply posts. View Active Threads