Reworking FFT code

[knc1]

Running the un-modified code under the DIY-KeK posted Aboriginal Linux using the DIY-KeK qemu-1.1.1 build; limited to a single cpucore - - -

Build notes:

(armv6l:1) fft-arm-0.01 # export CFLAGS="-O2 -march=armv6 -mfpu=vfp -mfloat-abi=softfp -fomit-frame-pointer"
(armv6l:1) fft-arm-0.01 # gcc $CFLAGS -o ffttest-arm radix4fft.c testfft.c testmain.c -lm
(armv6l:1) fft-arm-0.01 # strip --strip-unneeded ffttest-arm
(armv6l:1) fft-arm-0.01 # ls -l fft*
-rw-------    1 guest    guest         1959 Apr 17  2004 fft.h
-rwxr-xr-x    1 root     root         14068 Jul 17 12:08 ffttest-arm

Yields a test run of:

(armv6l:1) fft-arm-0.01 # ./ffttest-arm
Testing a 64-point FFT.
1:1 mapping; no reorder required.
SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 
5.87 us per 64-point FFT, or 21.80 Mbps with QPSK, or 20.29 insns per point on a 220MHz SA-1100.

NOTE: That hard coded machine type lies, this is an emulated armv6l running at 300MHz.

Seems like we could run one of those, with the most recent set of samples taken from the audio stream buffer, just before each e-ink update. Even at a blazing e-ink update rate of 15fps.

I.E: Adding another 6 us in front of each e-ink update to display a rt, graphic, equalizer with a spectral display isn't going to be noticeable - not even to the processor.

Let us see what the resource usage is when run as a stand-alone program (rather than hard coded into the equalizer):

(armv6l:1) fft-arm-0.01 # time ./ffttest-arm
Testing a 64-point FFT.
1:1 mapping; no reorder required.
SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 
5.90 us per 64-point FFT, or 21.71 Mbps with QPSK, or 20.38 insns per point on a 220MHz SA-1100.
real    0m 0.76s
user    0m 0.75s
sys    0m 0.01s

Hmm...
1 fps isn't going to keep up with a 15 (or 7) fps display.

Maybe hard coding the look-up tables ...

[knc1]

Originally Posted by knc1

Maybe hard coding the look-up tables ...

A first step in that direction, get a look at what the compiler is seeing:

(armv6l:1) fft-arm-0.01 # export CFLAGS="-E"
(armv6l:1) fft-arm-0.01 # gcc $CFLAGS -o ffttest-arm.i radix4fft.c testfft.c testmain.c

Skipping past all the "include *.h" bull ...
Now it can be seen (and read) just how those thousands of trig. calls are used to init the lookup tables.

The ffttest-arm.i compiler input file is gzip'd and attached here.

Oops
The -E option doesn't do multiple *.c files into a single listing.
Ah, but the reader should get the idea, even if it has to be called 3 times.

[geekmaster]

Originally Posted by knc1

...
Seems like we could run one of those, with the most recent set of samples taken from the audio stream buffer, just before each e-ink update. Even at a blazing e-ink update rate of 15fps.

I.E: Adding another 6 us in front of each e-ink update to display a rt, graphic, equalizer with a spectral display isn't going to be noticeable - not even to the processor.

The K3 does not return from the eink update calls until it is (usually) ready for you to write to the framebuffer again. How long it takes depends on the complexity of the scene change. For dithered animation the average speed is 7.7FPS (if you add a slight delay for less artifacts). If you go "by the book" you are supposed to wait longer (about 4 FPS) to avoid visible artifacts, and if you go even a little faster (7.8FPS) the animation quality does deteriorate significantly (which is why I settle on the 7.7FPS compromise). Although a K4/K5 can go faster with its different (SoC integrated) eink controller, it looks okay at 7.7FPS so I used that as my "standard" framerate in the GMV files.

The K4/K5 can run up to 45FPS with VERY carefully chosen animation, but for general unfiltered animation, 7.7FPS (average) is about as good as it gets. I can push out a little more speed without significant artifacts by doing spatiotemporal smoothing (averaging across multiple frames) to the video, which I did in a couple of tests, but that does not speed things up on the K3 so I do not make a habit of encoding with that smoothing.

It may turn out that how you display the FFT data may be able to do 15FPS on a K5 (and a K4 in diags mode), but the K4 in main mode emulate K3 eink (mostly) so it limits what kindle model you can use.

There is not much CPU left over on a K3 either while playing video with sound.

[knc1]

The one (working) FFT example from yesterday was written when "ARM doesn't have hardware floating point".

Our Kindles do have hardware floating point, almost, sort-of, ...
The K3 (VFP) unit **almost** supports IEEE 754 except for a few of the details that don't matter in the simple math world used here.

While re-vamping that decade old code (whose big claim to fame at the moment is that it "works out of the box") a person would want to keep in mind at least the basics:
http://infocenter.arm.com/help/index.../Chdidbba.html
(and the topic just above on the left)
That the VFP does best with short vectors of 8 single or 4 double precision values.

Also note, because of the **almost** IEEE 754, you want to include in the gcc optons:
-mfpu=vfp -mfloat-abi=softfp -funsafe-math-optimizations -std=c99
Also that libm is a double precision library (the -std=c99 is for some libm functions).

[knc1]

Now that we have a bit of a feel for how the emulated armv6l environment posted in DIY-KeK works (6us per 64 point FFT vs 9us per 64 point FFT on the K3) . . . .

Time to drop this gpl licensed fun show and go back to the thread topic, the bsd licensed one that needs to port from 16bit to 32 bit.

The 16bit example uses a DIY, 8Q7, real number storage format to fit the data into 16bits.
With our 32bit integer and floating point processors, we can us "int" and "float" (both 32bit register sizes) in place of the 8Q7.

Since libm does double precision (64bit) operations, we need to be a bit careful to keep local "float" versions of the libm results so that we don't need for the C code to "promote / demote" the data types a zillion times per math statement.

Although the VFP does have the float <-> double conversions in hardware, it is still better to avoid any work that can be avoided.

The VFP has 4 banks of (different views, same registers) either 8 single or 4 double precision vectors.
So we want to organize our progress as inner loops of either 4 or 8 data values (depending on precision at that place in the code).

And just so everyone does not have to go searching for the O.P. files, I have attached copies here:

[twobob]

This will clearly be awesome

[knc1]

A 64 band equalizer might well be considered the "minimum" for a sound studio.

When I had a nice listening room, I had a 16 band equalizer (with spectrum display - all electronic, not a computer anywhere).

Maybe setting our sights on an 8 band equalizer for the Kindles would be reasonable (after all, they are a "book" not a "studio").
People might expect to read about FFTs on a Kindle, not run one.

Plus, our hardware does math on 8 data point vectors per instruction.
Seems like that would fit with the overall plan.

It will probably take some clever coding of the overall algorithm to coax gcc into generating the code we want. That will be a learning experience, I am sure.

Let me see now,
we have 3, 8 data point vectors plus 8 float or integer registers for common constants.
So the question in my mind arises: "Can we do an 8 point FFT entirely in the VFP?"
Another interesting question to learn the answer to.

Hmm... I wonder if ARM has already published one as an "application note".
That one was too easy to answer:
http://infocenter.arm.com/help/index.../Cacjgfad.html
Now it just remains to find a copy - outside of the for-purchase development system that we can use. Or at least study for reference.

[twobob]

I have been digging around on this one.

http://hax4.blogspot.co.uk/2010_02_01_archive.html Advises:

"ARM has provided a folder under RVDS 4.0 installation path that contains the assembly code for FFT, which covers almost all the FFT scenarios, including the odd ones like radix = 3 (Can radix be 3 for FFT?). It supports both V4 and V5TE instruction set.

However it would be useful only if you know how to set it up. For some reason, ARM did not do a good job to document this assembly code. And it took me a while to sort it out :

(1) First of all, you need to determine the FFT length and the data format you are gonna use. Assume you are gonna use 16 bit complex number and N = 1024

For 16 bit complex number, each point should be like:

typedef struct { INT16 x, y; } comp16;


(2) Goto "C:\Program Files\ARM\RVDS\Examples\4.0\77\windows\fft_v5te", open the file fft_main.s under "fft_asm" folder, at the end of the file, you will find a bunch of "GENERATE_FFT_FUNCTION", which are macros to generate correspondent FFT functions for different scenarios. Only keep the one that fits your needs and comment out the rest. If will save you significant amount of memory space by doing do.

(3) Now you need to determine the function name for your FFT scenario. The assembly would generate your FFT function name in the following naming convention:

FFT_$radix.$direction.$qname.$intype.$outpos.$pres cale

For example, if your FFT is radix = 4, forward, 16 bit integer, different in/out buffer, no pre-scaling, your FFT function would be "FFT_4EFQ14". Please refer to fft_mac.h under "fft_asm" folder for detail of the parameter definition.

So your FFT function would be something like:

typedef int fFFT16(comp16 *input, comp16 *output, int N);
extern fFFT16 FFT_4EFQ14;
......
FFT_4EFQ14 (input_buffer, output_buffer, 1024); // actual FFT operation

(4) Also you need to generate the twiddle table for your FFT. Use gcc or VC++ to open the file gentable.c under "tables" folder, compile and run it. It will prompt you for the parameters regarding your FFT scenario, and generate a file like t_01024_4.c for you (depending on the parameters, the file name might be different). If you have a FFT with odd radix or length, you might also need to have the "t_rad.c" file.

Open the twiddle table you've just generated, only keep the table you that fits your FFT scenario and remove the rest. The table bears the name like t_QxxS_radix or t_QxxR_radix. Here the S and R have something to do with "coeforder", whose definition can be found in fft_mac.h as well. For our 1024 point, 16 bit FFT case, if the ARM architecture is no less than version 5, it should be R. So we should keep the table "t_Q14R_4" only.


Now put all the files together (all the files under "fft_asm" folder, plus the twiddle table you've generated) and you are good to go!

There is no doubt that those who composed this piece of FFT code are artful engineers. However, I do have a few notes on how the test program was written. After all the convoluted test set up in the testfft.c provided by ARM, I don't see how these "enum" type would help in any meaningful way. Correct me if I am wrong, but my understanding is that those "enum" types are mostly for the test program itself. Sadly the example test only gives out tests for floating format. I couldn't find any example that shows how to use 16 bit fixed-point.

Also the way function name is provided is kind sneaky. Although it is very clever to do this in assembly, such practice could really use some good documents."

I will follow up on this.
Thanks

[geekmaster]

Originally Posted by twobob

I have been digging around on this one.

http://hax4.blogspot.co.uk/2010_02_01_archive.html Advises:

I will follow up on this.
Thanks

Thanks for posting this. Regarding the radix-3 FFT, it was reported in some FFT code comments that an odd-radix FFT needs to be the last stage in a power-of-2 radix FFT pipeline. See one of the integer examples in the discussed FFT codees for the example of the 80-step FFT that has a radix-5 FFT as the last step.

[knc1]

Originally Posted by twobob

I have been digging around on this one.

http://hax4.blogspot.co.uk/2010_02_01_archive.html Advises:


I will follow up on this.
Thanks

Reads like it needs a strong dose of: gcc -E -o something.i ....
applied to it to turn it into readable "C" after all the set-up nonsense is done.

[knc1]

Here is a simple bit of math, good for someone (like me) who has never written VFP code.

We know that an array can be accessed by its member index.
We know that an array in C is just a pointer to the first (0) member.
We know that other fft implementations use either a lookup table or hard code something to "bit reverse" the member index numbers.
We know that C allows us to do "pointer math".
We know the VFP does not do "indexed load/store".

But it certainly does scaler, vector math.

We know the address of a member is: Maddr = (base + (index * (size_of(member)))
So how do we do "pointer math" for an 8 member array as a vector?

For an 8 element array, the bit reverse relationship looks like:
Fixed now
http://drpbox.knetconnect.com/fft/bits.html

Hmm (the intended algorithm, not the actual code) ...
* fill multiple with (size_of(member)) - bank 0
* load multiple with "magic index" - bank 1
* mul bank-0 bank-1 bank-0
* fill multiple with (array base) - bank 1
* add bank-0 bank-1 bank-0
and like magic we now have the pointer addresses in S0 .. S7 of the reordered input array.

C "pointer math" on steroids

[geekmaster]

Originally Posted by knc1

Here is a simple bit of math, good for someone (like me) who has never written VFP code.
...
C "pointer math" on steroids

That "broken image" in the post above is shown correctly in this embedded quote. Click this image to see the original HTML version.

[twobob]

Originally Posted by knc1

Here is a simple bit of math, good for someone (like me) who has never written VFP code.

We know that an array can be accessed by its member index.
We know that an array in C is just a pointer to the first (0) member.
We know that other fft implementations use either a lookup table or hard code something to "bit reverse" the member index numbers.
We know that C allows us to do "pointer math".
We know the VFP does not do "indexed load/store".

But it certainly does scaler, vector math.

We know the address of a member is: Maddr = (base + (index * (size_of(member)))
So how do we do "pointer math" for an 8 member array as a vector?

For an 8 element array, the bit reverse relationship looks like:
Fixed now
http://drpbox.knetconnect.com/fft/bits.html

Hmm (the intended algorithm, not the actual code) ...
* fill multiple with (size_of(member)) - bank 0
* load multiple with "magic index" - bank 1
* mul bank-0 bank-1 bank-0
* fill multiple with (array base) - bank 1
* add bank-0 bank-1 bank-0
and like magic we now have the pointer addresses in S0 .. S7 of the reordered input array.

C "pointer math" on steroids

Ooh. pretty...

It would seem my wine installation is Dead beyond recovery after several hours of digging around in registries. I have reported the bug (and found 4 other users who shared my experience, what a joy it is to belong)

To that end I shall be now wasting 4 hours doing this all again under Linux ; )

[knc1]

Originally Posted by geekmaster

That "broken image" in the post above is shown correctly in this embedded quote. Click this image to see the original HTML version.

The image embedded in your quote shows that I can not do add/subtract in my head.
The original is now corrected.

At least at re-vision #3 it is fixed, and now displays the symmetry I knew was in there.

[geekmaster]

Originally Posted by knc1

The image embedded in your quote shows that I can not do add/subtract in my head.
The original is now corrected.

I noticed that. The new text does not resize like the previous version that I "screenshotted". An HTML page does not work very well inside the IMG tags you have in your post.