Re: A "slanted edge" analysis program


Here is an implementation of Glassman's method of FFT, which wil work
for any N, not just powers of two. If N is not a power of two, it
degenerates to a DFT. The code is lifted from PAL (Public Ada
Library), if I remember correctly, and I do not think there is any
restrictions on it. You will have to convert it to C of course, but
the code is pretty obvious even if yo don't know Ada. Just beware that
unlike C, Ada allows nested subroutines, and that arrays do not
necessarily start with index 0.....

with Ada.Numerics.Short_Complex_Types;
use Ada.Numerics.Short_Complex_Types;

package FFT_Pack is

type Complex_Vector is array(Integer range <>) of Complex;

procedure Debug(X : Complex_Vector);

procedure FFT (FFT_Data : in out Complex_Vector ;
Inverse_Transform : in boolean );

end FFT_Pack;

with Ada.Numerics.Short_Elementary_Functions;
use Ada.Numerics.Short_Elementary_Functions;
with Ada.Text_Io; use Ada.Text_Io;

package body FFT_Pack is

procedure Debug(X : Complex_Vector) is
begin
for I in x'Range loop
Put_Line(Integer'Image(I) & " : "
& Short_Float'Image(X(I).Re) & " "
& Short_Float'Image(X(I).Im));
end loop;
end Debug;

procedure Glassman (A, B, C : in integer;
Data_Vector : in out Complex_Vector ;
Inverse_Transform : in boolean ) is

Temp : Complex_Vector(1..Data_Vector'length);
Counter : integer := Data_Vector'first;
JC : integer := 0;
Two_Pi : constant short_float := 6.28318530717958;
Del, Omega, Sum : Complex;
Angle : short_float;
C_Plus_1 : integer := C + 1;
begin
Temp(Temp'Range) := Data_Vector;
Angle := Two_Pi / (short_float(A * C));
Del := (Cos(Angle), (-(Sin(Angle))));
if (Inverse_Transform) then
Del := Conjugate(Del);
end if;
Omega := (1.0,0.0);

for IC in 1..C loop
for IA in 1..A loop
for IB in 1..B loop
Sum := Temp((((IA - 1)*C + (C-1))*B) + IB);
for JCR in 2..C loop
JC := C_Plus_1 - JCR; -- No need to add C + 1 each
-- time through loop
Sum := Temp((((IA - 1) * C + (JC - 1)) * B)
+ IB) + (Omega * Sum);
end loop; -- JCR
Data_Vector(Counter) := Sum;
Counter := Counter + 1;
end loop; -- IB
Omega := Del * Omega;
end loop; -- IA
end loop; -- IC
end Glassman;

procedure FFT ( FFT_Data : in out Complex_Vector;
Inverse_Transform : in boolean ) is

A : integer := 1;
B : integer := FFT_Data'length;
C : integer := 1;

begin -- FFT

while (B > 1) loop -- define the integers A, B, and C
A := C * A; -- such that A * B * C = FFT_Data'length
C := 2;
while (B mod C) /= 0 loop
C := C + 1;
end loop;
B := B/C; -- B = 1 causes exit from while loop
Glassman (A,B,C, FFT_Data, Inverse_Transform);
end loop;

if Inverse_Transform then -- optional 1/N scaling for inverse
-- transform only
for i in FFT_Data'range loop
FFT_Data(i) := FFT_Data(i) /
short_float(FFT_Data'length);
end loop;
end if;
end FFT;
end FFT_Pack;


--
C++: The power, elegance and simplicity of a hand grenade.
.

Relevant Pages

  • Re: Idle Curiosity -- 4-20mA Signalling
    ... 4-20mA signaling loop? ... I'm curious as to how much freedom one has to power ones device when one ... designs some gizmo that flaps in the breeze on the end of a current loop. ... maximum supply voltage and calculate the maximum load resistance. ...
    (sci.electronics.design)
  • Re: FFT for N Not a Power of 2
    ... remove the largest number of observations that are a power of 2, ... There are FFT algorithms that work for any N, including prime numbers, but ... with floor rising near the two tones. ... had to ensure that the length of the window and periods of each of the ...
    (sci.optics)
  • Re: More on cheater plugs
    ... Two hots for POWER and a good ground. ... interconnects, that have their 3rd wire safety pins connected to each other, ... BUT, for the sake of this gedanken experiment, the 3rd wire safeties are not ... There is a closed loop path for electrons. ...
    (rec.audio.opinion)
  • Re: Battery Charge Noise/Interference on Audio Output
    ... It is in fact noise caused by an 'earth loop'. ... from the power supply. ... This depends on where you experienced your earth loop. ...
    (comp.sys.laptops)
  • Re: Facing Problem in Loop
    ... it goes to infinite loop. ... Second problem is it works fine for smaller value of X or N but not ... for higher say base 15 power 20. ... double base, pwr; ...
    (comp.lang.c)