I’m doing some frequency analysis - as usual I trying to minimise the processing cycles to minimise current consumption.
I’ve only just discovered the Fast Hartley Transform (FHT). It seems pretty good ( and simple) - I wonder why it is not more widely used? Anyway, maybe others can give it a go - here is the info to get started:
I got the Pascal version of the FHT from the following document-
“OPTIMIZED FAST HARTLEY TRANSFORM FOR THE MC68000 WITH APPLICATIONS IN IMAGE PROCESSING A THESIS SUBMITTED TO THE FACULTY IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE BY
A. ARLO REEVES�?
On page 93 there is a Pascal program under the following title-
“DFHT - DFHT is a high level implementation of the FHT, included here for those who want more than just assembly source. .....�?
The program is easily converted to C , and the routine is compact enough to convert into a single assembly language routine . My first attempt with the Cortex-M4 (using F32 with the FPU) resulted in a routine that uses 104, 000 cycles for a 1024 point FHT- much faster than my attempts using a radix2 FFT. Moving to fixed point arithmetic later will give even more savings.
However for testing the above FHT program, the following two routines are needed - one for converting the output of the FHT into amplitudes, and the other routine for converting the output of the FHT to the same output as from a normal FFT. You can use the output data from these routines (using pseudo random input) to verify your program against the output of an online FFT tool. I found the following one very useful - http://scistatcalc.blogspot.com/2013/12/fft-calculator.html
procedure map_to_fourier;
var
i, n_minus_i : integer;
begin
for i := 0 to maxN - 1 do
begin
if i = 0 then
n_minus_i := 0
else
n_minus_i := maxN - i;
reals [ i ] := ( x[ i ] + x[ n_minus_i ]) / 2;
imags [ i ] := - ( x[ i ] - x[ n_minus_i ]) / 2;
end;
end;
procedure calc_power;
var
i, n_minus_i : integer;
begin
for i := 0 to maxN - 1 do
pwr[i] := 0;
for i := 0 to (maxN div 2) do
begin
if i = 0 then
n_minus_i := 0
else
n_minus_i := maxN - i;
pwr [i] := sqrt( (x[i]*x[i] + x[n_minus_i]*x[n_minus_i]) / 2 );
end;
end;
