RFC: Generic Fixed-Point IIR Filter

RFC: Generic Fixed-Point IIR Filter

[Patrick Noffke]

I would greatly appreciate some feedback/help on implementing a generic fixed-point IIR filter package.

I am attempting to define a fixed-point type (Pressure_Type) that has a limited range and precision.  When doing things (like filtering) to values of this type, intermediate calculations require a larger range and possibly higher precision.  But I don't want to expand the range/precision of the base type -- rather define new types as appropriate for the required intermediate calculations.

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?

 - 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?

 - I get warnings that the filter constants are not a multiple of small.  I understand why this is (small is the largest (most positive) negative power of two not greater than delta).  I can explicitly modify the constants to be a multiple of small, but is there another way to define the constants that will "convert" them to a multiple of small?

 - I realize my range and precision for Filter_Value_Type are not necessarily ideal (the range depends on the filter coefficients, the filter architecture, as well as the filter order).  Do you know of any references or have suggestions for how to analyze a filter to determine the ideal range and precision?


iir_filter.ads:

generic
   --  The fixed-point numeric type for the filter input and output.
   type Fixed_Type is delta <>;
   --  The fixed-point numeric type for intermediate filter calculations.
   --  Should have enough range to avoid overflow and enough precision so
   --  rounding errors are insignificant.
   type Filter_Value_Type is delta <>;
   --  Number of "a" coefficients.
   N_A : Positive;
   --  Number of "b" coefficients.
   N_B : Positive;
package IIR_Filter is
   type A_Coefficient_Array_Type is array (1 .. N_A) of Filter_Value_Type;
   type B_Coefficient_Array_Type is array (1 .. N_B) of Filter_Value_Type;
   type IIR_Filter_Type is private;

   procedure Init_Filter (Filter         : in out IIR_Filter_Type;
                          A_Coefficients :        A_Coefficient_Array_Type;
                          B_Coefficients :        B_Coefficient_Array_Type;
                          Initial_Value  :        Fixed_Type);

   function Step_Filter (Filter : in out IIR_Filter_Type;
                         X      :        Fixed_Type) return Fixed_Type;

private
   type Y_Array_Type is array (1 .. N_A) of Filter_Value_Type;
   type X_Array_Type is array (1 .. N_B) of Filter_Value_Type;
   type IIR_Filter_Type is record
      A : A_Coefficient_Array_Type;
      B : B_Coefficient_Array_Type;
      X_Array : X_Array_Type;
      Y_Array : Y_Array_Type;
   end record;

end IIR_Filter;

iif_filter.adb:

package body IIR_Filter is

   procedure Init_Filter (Filter         : in out IIR_Filter_Type;
                          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;
      for Y_Index in Filter.Y_Array'Range loop
         Filter.Y_Array (Y_Index) := Filter_Value_Type (Initial_Value);
      end loop;
   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
         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;

      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;

      --  Compute -a(2) * y(n - 1) - a(3) * y(n - 2) + ...
      if Filter.Y_Array'Length > 1 then
         for Y_Index in 2 .. Filter.Y_Array'Last loop
            Filter.Y_Array (1) := Filter.Y_Array (1) -
              Filter.A (Y_Index) * Filter.Y_Array (Y_Index);
         end loop;
      end if;

      --  Divide by a(1).
      Filter.Y_Array (1) := Filter.Y_Array (1) / Filter.A (1);

      return Fixed_Type (Filter.Y_Array (1));
   end Step_Filter;

end IIR_Filter;

pressure.ads:

package Pressure is
   type Pressure_Type is delta 0.001 range -20.0 .. 400.0;
end Pressure;


Then in my main procedure:

procedure Filter_Test is
   N_A : constant Positive := 3;
   N_B : constant Positive := 3;
   type Filter_Value_Type is delta Pressure.Pressure_Type'Delta / 10
     range N_A * Pressure.Pressure_Type'First .. N_A * Pressure.Pressure_Type'Last;
   package Pressure_IIR_Filter is new
     IIR_Filter (Fixed_Type => Pressure.Pressure_Type, Filter_Value_Type => Filter_Value_Type,
                 N_A => N_A, N_B => N_B);
   Pressure_Filter : Pressure_IIR_Filter.IIR_Filter_Type;
   A_Coeffs : constant Pressure_IIR_Filter.A_Coefficient_Array_Type :=
     (1.0, -1.84036, 0.85219);
   B_Coeffs : constant Pressure_IIR_Filter.B_Coefficient_Array_Type :=
     (0.0029583, 0.0059165, 0.0029583);
   Filtered_Pressure : Pressure.Pressure_Type := 0.0;

begin

   Pressure_IIR_Filter.Init_Filter (Pressure_Filter, A_Coeffs, B_Coeffs, 0.0);

   for Step_Index in 1 .. 200 loop
      Filtered_Pressure :=
        Pressure_IIR_Filter.Step_Filter (Pressure_Filter, 300.0);
      Put_Line ("Step: " & Step_Index'Img & ", Filtered_Pressure = " &
                  Filtered_Pressure'Img);
   end loop;

end Filter_Test;


Regards,
Patrick

Re: RFC: Generic Fixed-Point IIR Filter

[Dmitry A. Kazakov]

On Thu, 26 Mar 2015 09:25:25 -0700 (PDT), Patrick Noffke wrote:

> I realize my range and precision for Filter_Value_Type are not
> necessarily ideal (the range depends on the filter coefficients, the
> filter architecture, as well as the filter order).  Do you know of any
> references or have suggestions for how to analyze a filter to determine
> the ideal range and precision?

Intermediate values should have some physical sense. E.g. acceleration when
you calculate dv/dt from v and t. Such values are limited by the process
and thus you could estimate their ranges.

This is also a good argument against using a single fixed-point type in
calculations. This is not a problem with floating-point types as they
designed to accommodate all possible intermediate values of practically any
process, though maybe losing some accuracy. But with fixed-point types one
should go after proper dimensional arithmetic with *different* types of
*different* precision for each physical value, intermediates included.

Unfortunately there is basically no support in Ada (or in any other
language I know) for this.

-- 
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de

Re: RFC: Generic Fixed-Point IIR Filter

[Jeffrey Carter]

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

Re: RFC: Generic Fixed-Point IIR Filter

[Patrick Noffke]

On Thursday, March 26, 2015 at 12:41:29 PM UTC-5, Jeffrey Carter wrote:
> On 03/26/2015 09:25 AM, Patrick Noffke wrote:
> 
> 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.
> 

I didn't know you could do this.  Thanks.

> > 
> >  - 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 should have known! :-)

> 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,

How is universal_fixed (likely) implemented on a 32-bit processor?  What happens if you multiply two fixed-point types that use all 32-bits for range and precision?  What if one type uses all the bits for the range and another uses all the bits for precision?

I'm targeting an ARM Cortex-M4 processor, but writing test code on an x86_64 PC.  In my test code with FSF GNAT, I noticed from a hex print of a Pressure_Type variable that only the minimum number of bits were used (as I understand it, with Ada 95 and later, small *could* be smaller than the maximum allowable small).

I know the M4 processor has the SMULL instruction (Signed Multiply (32x32), 64-bit result) and similar instructions, but I don't yet know if the compiler is using these for 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.
> 

True, for the generic package.  But I really only need a 2nd order filter for my application, so I should be able to do:

Filter.Y_Array (1) := (Filter.B (1) * X (1) + Filter.B (2) * X (2) + Filter.B (3) * X (3) - Filter.A (2) * Filter.Y_Array (2) - Filter.A (3) * Filter.Y_Array (3)) / Filter.A (1);

Then every product is universal_fixed, so the sum is universal_fixed, and implicitly converted to the LHS type.  But again, I'd like to know what is happening under the hood to implement this.

> Some comments on your implementation:
> 

Thank you for those.  I will address the issues you point out.

Patrick

Re: RFC: Generic Fixed-Point IIR Filter

[Randy Brukardt]

"Patrick Noffke" <patrick.noffke@gmail.com> wrote in message 
news:2c2ee86e-b9bd-49e3-aa7f-206f3c4da95e@googlegroups.com...
...
> How is universal_fixed (likely) implemented on a 32-bit processor?
> What happens if you multiply two fixed-point types that use all 32-bits
> for range and precision?  What if one type uses all the bits for the range
> and another uses all the bits for precision?

There's probably not a single implementation; the implementation likely 
depends on the smalls and ranges of the input types and the output 
(converting) type. (There has to be a convert-to type, universal fixed can 
never stand alone.) In particular, in simple cases (where there is a clear 
relationship between the small values such that scaling is an integer 
operation), probably a 32-bit or 64-bit integer operation would be used. In 
cases that aren't so simple, who knows? Janus/Ada falls back on floating 
point in that case, because what do you do if the small of one input type is 
0.3141596 and the other is 0.0025, and the small of the output is 0.0001? 
Hopefully, such things don't occur in code written by people that know what 
they're doing, but it's possible and the compiler has to somehow generate 
code for it. At least the floating point temporaries have lots of bits 
available, so you're probably going to get the right answer if possible.

Anyway, you're better off reasoning from principles. Universal fixed is not 
supposed to overflow (it doesn't have a constraint), so I would expect some 
sort of 64-bit implementation ultimately.

Note that this is wrong:

>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.

A result of universal-fixed always has to be converted immediately to 
another type. It's not legal to use it directly. Specifically, 4.5.5(19.1/2) 
says: "The above two fixed-fixed multiplying operators shall not be used in 
a context where the expected type for the result is itself universal_fixed 
".

That's so the compiler can always know how many bits are needed for the 
intermediate results (in order to avoid overflow). But you'll need to figure 
those out yourself (and as you noted, you can't declare them inside of the 
generic).

There's a reason that fixed point is rarely used!

                                        Randy.



                                     Randy.

Re: RFC: Generic Fixed-Point IIR Filter

[Jeffrey Carter]

On 03/26/2015 12:44 PM, Patrick Noffke wrote:
> 
> How is universal_fixed (likely) implemented on a 32-bit processor?  What
> happens if you multiply two fixed-point types that use all 32-bits for range
> and precision?  What if one type uses all the bits for the range and another
> uses all the bits for precision?

I refer you to ARM 3.4.1(7):

"The set of values of a universal type is the undiscriminated union of the set
of values possible for any definable type in the associated class." So
universal_fixed has to be able to handle any value of any fixed-point type
definition.

-- 
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

Re: RFC: Generic Fixed-Point IIR Filter

[Jeffrey Carter]

On 03/26/2015 01:34 PM, Randy Brukardt wrote:
> 
> A result of universal-fixed always has to be converted immediately to 
> another type. It's not legal to use it directly. Specifically, 4.5.5(19.1/2) 
> says: "The above two fixed-fixed multiplying operators shall not be used in 
> a context where the expected type for the result is itself universal_fixed 
> ".

So there's always an expected subtype? I thought in a string of operations in an
expression, there was no expected type until the conversion of the result to the
expected subtype of the entire expression.

-- 
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

Re: RFC: Generic Fixed-Point IIR Filter

[Randy Brukardt]

"Jeffrey Carter" <spam.jrcarter.not@spam.not.acm.org> wrote in message 
news:mf1qo9$dqb$3@dont-email.me...
> On 03/26/2015 01:34 PM, Randy Brukardt wrote:
>>
>> A result of universal-fixed always has to be converted immediately to
>> another type. It's not legal to use it directly. Specifically, 
>> 4.5.5(19.1/2)
>> says: "The above two fixed-fixed multiplying operators shall not be used 
>> in
>> a context where the expected type for the result is itself 
>> universal_fixed
>> ".
>
> So there's always an expected subtype? I thought in a string of operations 
> in an
> expression, there was no expected type until the conversion of the result 
> to the
> expected subtype of the entire expression.

Right. The rule specifically does not allow something like:
    A * B * C
because that could take an infinite number of bits to evaluate (just keeping 
adding the * operators) - the number of bits needed is the sum of the 
maximum bits for each value. Universal-fixed is not supposed to overflow, 
some something has to make it implementable. In addition, most machines have 
a double precision multiply operation (like the 32*32 bits gives 64 bits 
operation the OP mentioned and a matching divide). So there's an obvious 
implementation in many cases (not on the margins, as I noted), so long as we 
are talking about a single multiply. Ada runtime math isn't supposed to 
require some complex infinite precision operation (which admittedly is 
interesting, since the compiler is required to do exactly that).

                                     Randy.

Re: RFC: Generic Fixed-Point IIR Filter

[Patrick Noffke]

On Thursday, March 26, 2015 at 3:34:23 PM UTC-5, Randy Brukardt wrote:
> "Patrick Noffke" wrote in message 
> ...
> > How is universal_fixed (likely) implemented on a 32-bit processor?
> > What happens if you multiply two fixed-point types that use all 32-bits
> > for range and precision?  What if one type uses all the bits for the range
> > and another uses all the bits for precision?
> 
> Note that this is wrong:
> 
> >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.
> 
> A result of universal-fixed always has to be converted immediately to 
> another type. It's not legal to use it directly. Specifically, 4.5.5(19.1/2) 
> says: "The above two fixed-fixed multiplying operators shall not be used in 
> a context where the expected type for the result is itself universal_fixed 
> ".
> 

When I try to compute a 2nd order filter result as:

            Filter.Y (1) := (Filter.B (1) * Filter.X (1) +
                               Filter.B (2) * Filter.X (2) +
                               Filter.B (3) * Filter.X (3) -
                               Filter.A (2) * Filter.Y (2) -
                               Filter.A (3) * Filter.Y (3)) / Filter.A (1);

I get this error:

iir_filter.adb:46:60: ambiguous universal_fixed_expression
iir_filter.adb:46:60: possible interpretation as type "Standard.Duration"
iir_filter.adb:46:60: possible interpretation as type "Fixed_Type" defined at iir_filter.ads:3


I don't really understand the clause in 4.5.5(19.1/2).  What defines the context?  I think it's clear in my statement the "expected type for the result" is Fixed_Type, so why is Jeffrey's statement wrong?

And then how do I get past this error?  Is the problem with the + operator?  How do I remove the ambiguity?

> That's so the compiler can always know how many bits are needed for the 
> intermediate results (in order to avoid overflow). But you'll need to figure 
> those out yourself (and as you noted, you can't declare them inside of the 
> generic).
> 
> There's a reason that fixed point is rarely used!
> 

My processor doesn't have an FPU.

Patrick

Re: RFC: Generic Fixed-Point IIR Filter

[Patrick Noffke]

On Thursday, March 26, 2015 at 4:12:30 PM UTC-5, Randy Brukardt wrote:
> Right. The rule specifically does not allow something like:
>     A * B * C

Ah, okay.

If I change the calculation to

            Filter.Y (1) := Filter.B (1) * Filter.X (1) +
              Filter.B (2) * Filter.X (2) +
              Filter.B (3) * Filter.X (3) -
              Filter.A (2) * Filter.Y (2) -
              Filter.A (3) * Filter.Y (3);
            Filter.Y (1) := Filter.Y (1) / Filter.A (1);

Then it compiles.  So it didn't like me doing the division with a universal_fixed product as the numerator.

In my case, A(1) is 1.0, but that is not always the case, so there could still be overflow.  I suppose I could change this to a Normalized_IIR_Filter (implying A(1) is 1.0).

Thanks for your help.

Patrick

Re: RFC: Generic Fixed-Point IIR Filter

[Jeffrey Carter]

On 03/26/2015 02:12 PM, Randy Brukardt wrote:
> 
> Right. The rule specifically does not allow something like:
>     A * B * C
> because that could take an infinite number of bits to evaluate (just keeping 
> adding the * operators) - the number of bits needed is the sum of the 
> maximum bits for each value. Universal-fixed is not supposed to overflow, 
> some something has to make it implementable. In addition, most machines have 
> a double precision multiply operation (like the 32*32 bits gives 64 bits 
> operation the OP mentioned and a matching divide). So there's an obvious 
> implementation in many cases (not on the margins, as I noted), so long as we 
> are talking about a single multiply. Ada runtime math isn't supposed to 
> require some complex infinite precision operation (which admittedly is 
> interesting, since the compiler is required to do exactly that).

That makes sense. The main reason for fixed point in Ada 83 was to be faster
than SW floating-point on processors without FPUs. But since, as you note, the
compiler has to be able to do unbounded-precision math, then the compiler could
link in the library it uses to do that and do its intermediate calculations
using it. There could be a compiler switch to determine which is done. I think
GNAT has this.

What would be even better would be for the language to make the spec of that
unbounded-precision package part of the standard library.

-- 
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

Re: RFC: Generic Fixed-Point IIR Filter

[Brian Drummond]

On Thu, 26 Mar 2015 09:25:25 -0700, Patrick Noffke wrote:

> I would greatly appreciate some feedback/help on implementing a generic
> fixed-point IIR filter package.
...
> 
>  - 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?

I wrote (adapted) a generic SVD (matrix operation) to instantiate with 
either fixed or floating point. The original was floating point : the 
intent was to investigate if different fixed point formats could be made 
stable, before porting to VHDL and implementing in hardware. It quickly 
answered the question : for matrix inversion, no way!

Intermediate addition operations seemed OK : intermediate operations 
involving multiplication (or division) had to be rewritten using 
intermediate variables. 

A moment's thought showed this made sense : you may want the intermediate 
results stored at a different resolution : this is also usually the case 
with IIR filters.

In other words I had to break up 
   f := a * b * c;
into something like
   temp := b * c; 
   f := a * temp;
where temp may have been a higher resolution type, which allowed (i.e. 
forced!) me to explicitly control the numerical properties. This is the 
point Randy made, 4.5.5(19.1/2).

>  - is there another way to
>  define the constants that will "convert" them to a multiple of small?

I would probably use a conversion function on real (float) constants.

>  - I realize my range and precision for Filter_Value_Type are not
>  necessarily ideal (the range depends on the filter coefficients, the
>  filter architecture, as well as the filter order).  Do you know of any
>  references or have suggestions for how to analyze a filter to determine
>  the ideal range and precision?

Signal processing texts aimed at hardware engineers should cover this : 
"Synthesis and optimisation of DSP Algorithms" (Constantinides,Cheung,Luk) 
comes to mind. But with a generic, you can instantiate it with different 
fixed (and float) types, measure the quality of each, and choose.

-- Brian

Re: RFC: Generic Fixed-Point IIR Filter

[G.B.]

On 26.03.15 22:12, Randy Brukardt wrote:
> Ada runtime math isn't supposed to
> require some complex infinite precision operation (which admittedly is
> interesting, since the compiler is required to do exactly that).

The requirement to not do infinite precision arithmetic
for Standard."*" seems plausible to me insofar as Standard."*"
in source text conveys "a single multiplication", not
some big number operation which could take unpredictable
amounts of resources. (Still interesting that compiler
capabilities are withheld from its users, though.)