Re: RFC: Generic Fixed-Point IIR Filter

[Jeffrey Carter <spam.jrcarter.not@spam.not.acm.org>]

On 03/26/2015 09:25 AM, Patrick Noffke wrote:
> 
> Below is my first attempt.  I noticed the following issues:
>  - I cannot define Filter_Value_Type within the IIR_Filter package (i.e. it has to be a generic parameter).  If I try to put it in the package and base it on Fixed_Type, I get an error that a non-static expression is used for the delta and range.  For example, this gives the error I just mentioned:
> 
>    type Filter_Value_Type is delta Fixed_Type'Delta / 10
>      range N_A * Fixed_Type'First .. N_A * Fixed_Type'Last;
> 
> I realize (from http://en.wikibooks.org/wiki/Ada_Programming/Generics) that generic formal objects are non-static, but it's not clear to me why this won't work for a generic formal type (how can a type be non-static?).  Is there another way to parameterize the delta and range for intermediate calculations?

All generic formal parameters are non-static. As to a non-static type consider

generic -- P
   type T is range <>;
package P is ...

and then

function Input return Natural;
-- Obtains a value of subtype Natural from the user

...

subtype U is Integer range 1 .. Input;

package Q is new P (T => U);

Clearly subtype U is non-static, and so formal type T is also non-static.

> 
>  - If I try to use Fixed_Type for X_Array and Y_Array, then I do in fact get a constraint error (due to range overflow for an intermediate calculation).  Is there a better way to do a calculation for a fixed-point type when you expect the result to be within the specified range of Fixed_Type, but intermediate calculations may not?

All your problems result from suffixing type names with _Type :)

I refer you to ARM 4.5.5, specifically ¶18-19. Multiplication and division
between fixed-point values gives results of type universal_fixed, which is
implicitly convertible to any fixed-point type. When you write

Filter.Y_Array (1) := Filter.B (1) * X;

you are multiplying a Filter_Value_Type * Fixed_Type yielding universal_fixed,
which is then implicitly converted to the type of the LHS, Filter_Value_Type.
Your exceptions result from attempting to convert intermediate results to
Fixed_Type, which lacks the range needed to hold the values. It is therefore
sometimes possible to avoid the 2nd type if all calculations can be done in a
single expression which yields a value that fits in the smaller type, despite
having intermediate results which do not. That doesn't seem to be case for your
problem, though.

Some comments on your implementation:

> package body IIR_Filter is
> 
>    procedure Init_Filter (Filter         : in out IIR_Filter_Type;

Filter appears to be an "out" parameter, not "in out".

>                           A_Coefficients :        A_Coefficient_Array_Type;
>                           B_Coefficients :        B_Coefficient_Array_Type;
>                           Initial_Value  :        Fixed_Type) is
>    begin
>       Filter.A := A_Coefficients;
>       Filter.B := B_Coefficients;
>       for X_Index in Filter.X_Array'Range loop
>          Filter.X_Array (X_Index) := Filter_Value_Type (Initial_Value);
>       end loop;

Filter.X_Array := (others => Filter_Value_Type (Initial_Value) ); ?

>       for Y_Index in Filter.Y_Array'Range loop
>          Filter.Y_Array (Y_Index) := Filter_Value_Type (Initial_Value);
>       end loop;

Or even

Filter := (A       => A_Coefficients,
           B       => B_Coefficients,
           X_Array => (others => Filter_Value_Type (Initial_Value) ),
           Y_Array => (others => Filter_Value_Type (Initial_Value) ) );

>    end Init_Filter;
> 
>    function Step_Filter (Filter : in out IIR_Filter_Type;
>                          X      :        Fixed_Type) return Fixed_Type is
>    begin
>       --  Compute IIR filter using Direct Form I.
> 
>       --  Move the elements to next oldest position.
>       --  TODO:  Use a circular buffer for this.
>       if Filter.X_Array'Length > 1 then

This if statement (and others like it) is unnecessary. The loop will execute
zero times if the length is 1.

>          for X_Index in reverse 2 .. Filter.X_Array'Last loop
>             Filter.X_Array (X_Index) := Filter.X_Array (X_Index - 1);
>          end loop;
>       end if;

Filter.X_Array (2 .. Filter.X_Array'Last) :=
   Filter.X_Array (1 .. Filter.X_Array'Last - 1);

> 
>       if Filter.Y_Array'Length > 1 then
>          for Y_Index in reverse 2 .. Filter.Y_Array'Last loop
>             Filter.Y_Array (Y_Index) := Filter.Y_Array (Y_Index - 1);
>          end loop;
>       end if;
> 
>       -- Store the new X value.
>       Filter.X_Array (1) := Filter_Value_Type (X);
> 
>       --  Compute new filter output.  Initialize with b(1) * x(n).
>       Filter.Y_Array (1) := Filter.B (1) * X;
> 
>       --  Compute b(2) * x(n - 1) + b(3) * x(n - 2) + ...
>       if Filter.X_Array'Length > 1 then
>          for X_Index in 2 .. Filter.X_Array'Last loop
>             Filter.Y_Array (1) := Filter.Y_Array (1) +
>               Filter.B (X_Index) * Filter.X_Array (X_Index);
>          end loop;
>       end if;

You appear to be storing things backwards: x(n) is in X_Array (1), x(n-1) in
X_Array (2), ... (similarly for y). This makes your comments confusing.

HTH.

-- 
Jeff Carter
"Since I strongly believe that overpopulation is by
far the greatest problem in the world, this [Soylent
Green] would be my only message movie."
Charleton Heston
123