From: albinali on
Hi,

I am trying to implement the filter function in Matlab and it is
described as "The filter function is implemented as a direct form II
transposed structure:

y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
- a(2)*y(n-1) - ... - a(na+1)*y(n-na)

where n-1 is the filter order, which handles both FIR and IIR filters,
na is the feedback filter order, and nb is the feedforward filter
order.

I wrote the following code to implement it but I am still getting
different results. My filtered output seems to be monotonically
increasing. I would appreciate it if someone can take a look at my
code and let me know what I did wrong. Thanks!

Fahd

static double[] b = new double[5] {1.0,0.0,-2.0,0.0,1.0};
static double[] a = new double[5] { 1.00, -3.9946673,
5.9840174, -3.9840329, 0.9946827 };

static int NZEROS = 4;
static int NPOLES = 4;
static int ORDER = 4;
static double[] xv = new double[NZEROS + 1];
static double[] yv = new double[NPOLES + 1];
static double[] Filter(double[] data)
{

double[] filtered = new double[data.Length];

for (int i = 0; (i < data.Length); i++)
{
xv[0] = xv[1]; xv[1] = xv[2]; xv[2] = xv[3]; xv[3] =
xv[4];
xv[4] = data[i];
yv[0] = yv[1]; yv[1] = yv[2]; yv[2] = yv[3]; yv[3] =
yv[4];
yv[4] = b[0] * xv[4] + b[1] * xv[3] + b[2] * xv[2] +
b[3] * xv[1] + b[4] * xv[0] -
a[1] * yv[3] - a[2] * yv[2] - a[3] * yv[1] - a[4]
* yv[0];


filtered[i] = yv[4];
}

return filtered;
}
From: Tim Wescott on
On 05/25/2010 09:21 AM, albinali wrote:
> Hi,
>
> I am trying to implement the filter function in Matlab and it is
> described as "The filter function is implemented as a direct form II
> transposed structure:
>
> y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
> - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
>
> where n-1 is the filter order, which handles both FIR and IIR filters,
> na is the feedback filter order, and nb is the feedforward filter
> order.
>
> I wrote the following code to implement it but I am still getting
> different results. My filtered output seems to be monotonically
> increasing. I would appreciate it if someone can take a look at my
> code and let me know what I did wrong. Thanks!
>
> Fahd
>
> static double[] b = new double[5] {1.0,0.0,-2.0,0.0,1.0};
> static double[] a = new double[5] { 1.00, -3.9946673,
> 5.9840174, -3.9840329, 0.9946827 };

First, the unstable pole at z = 1.016 may have something to do with the
monotonic (and no doubt exponential) increase.

Second, as polynomial order goes up the roots become exceedingly
sensitive to variations in the coefficients. This creates all sorts of
numerical havoc with direct-form IIR filters. Thus, any time you
implement an IIR filter you really, _really_ want to factor it as far
down as you can, with second-order filters for the complex pole pairs,
with first-order filters for any real poles. Then implement your
complete filter as a cascade of the individual filters.

Second and a half, your unstable pole may be an unintended consequence
of this sensitivity, and the fact that you're expressing the filter as a
4-th order polynomial.

HTH

--
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com
From: Mikolaj on
Code looks fine.
Your filter is unstable
so your results are fine.