FIR Bandpass filtering (C code)

Hi folks,

I'm trying to implement an FIR bandpass filter in C; my reference material
is the excellent "Understanding Digital Signal Processing (2nd ed)" by
Richard G Lyons.

I'm starting with the lowpass filter, and then modifying the parameters for
bandpass (described below).

I'm calculating the lowpass filter coefficients based on the formula
on page 169:

        1    sin (PI * k * K / N)
h (k) = -- x --------------------
        N      sin (PI * k / N)

I'm filtering stock data, so my sampling frequency is one sample per day,
and I'm interested in bandpass filtering 8 day periods, 16 day, 32 day, and
higher, right up to multi-year trends.

I pass my function ("compute_filter()") the period (number of days) and
the filter bandwidth in days (period, rather than frequency).  For
example, "8" and "0.8", hoping for a 8 +/- 0.8 day filter...
The function returns the coefficients array and the number of elements.

The code calculates the lowpass formula above, and then multiplies the
coefficients by "a sinusoid of frequency X" (calculated as 1 / period)
as per "Bandpass FIR Filter Design" on page 183:

hbp (k) = hlp (k) * sshift (k)

Here's the code for reference -- I allocate the "f" array (the coefficients)
dynamically and size it dynamically.  This is done in an attempt to make the
coefficients array as small as required (the computation of N).

static double *
compute_filter (int *ntaps, int period, double bw)
{
    double      *f;         // the output array
    double      fm;
    double      x;
    double      N;
    int         nt2;        // N / 2 as an integer
    int         k;
    double      K;

    fm = 1. / (double) period;

    // ensure we have an odd number of samples
    N = (double) period * 2. + 1.;
    nt2 = N / 2;

    K = N / bw + 1.;

    f = calloc ((int) N, sizeof (*f));
    // calloc error checking deleted

    // calculate the lowpass coeff's
    for (k = 1; k <= nt2; k++) {
        x = M_PI * (double) k / N;
        f [nt2 + k] = f [nt2 - k] = sin (K * x) / sin (x) / N;
    }
    f [nt2] = K / N;

    // multiply by sinusoid to shift frequency center
    for (k = 1; k <= nt2; k++) {
        f [nt2 + k] *= cos (fm * (double) k * 2. * M_PI);
        f [nt2 + k] *= cos (fm * (double) k * 2. * M_PI);
    }
    // f [0] stays as is

    *ntaps = N;
    return (f);
}

The problem I'm running into is that the output of the filter follows the
input shape (the low frequency component).  Filtering seems to be happening,
in that the higher frequencies are being filtered out (the graph looks
smoother and smoother as the number of days is increased).  The filter
is acting more like a lowpass filter, rather than bandpass.

The convolution is done in another routine by:

    // n is the # of coeffs
    // f is the filter coeff array
    for (i = 0; i < ntrades - n; i++) {
        sum = 0;
        for (j = 0; j < n; j++) {
            sum += f [j] * closing_prices [i + j];
        }
        printf ("%d %g\n", i + n/2, sum);
    }

A few questions:
    a) how should the sin (k*K)/sin(k) function avoid a divide by zero?
       (I got around it by starting the array calculation at 1, and
       calculating the [0] element manually.)
    b) is my calculation of the minimum number of coefficients correct?
       Do I need 2X the number of coeffs in ALL cases, or can I get away
       with less in low frequencies or high frequencies?
    c) is my assumption about the bandwidth size correct (i.e., 8 day
       period with a 0.8 day bandwidth)?
    d) why doesn't it work? ;-)

Thank you very much in advance!

Cheers,
-RK
-- 
Robert Krten, Antique computer collector looking for PDP-8 and PDP-8/S
minicomputers; check out their "good home" at www.parse.com/~museum

On 1 Apr, 04:33, info2...@parse.com (Robert Krten) wrote:
> Hi folks,
>
> I'm trying to implement an FIR bandpass filter in C; my reference material
> is the excellent "Understanding Digital Signal Processing (2nd ed)" by
> Richard G Lyons.
>
> I'm starting with the lowpass filter, and then modifying the parameters for
> bandpass (described below).
>
> I'm calculating the lowpass filter coefficients based on the formula
> on page 169:
>
>         1    sin (PI * k * K / N)
> h (k) = -- x --------------------
>         N      sin (PI * k / N)
>
> I'm filtering stock data, so my sampling frequency is one sample per day,
> and I'm interested in bandpass filtering 8 day periods, 16 day, 32 day, and
> higher, right up to multi-year trends.

One cycle / year = 1/ 300 or whatever the number of trading
days per year is... narrow bands. You might need to get
into advanced techniques to get this one to work.

...


> The problem I'm running into is that the output of the filter follows the
> input shape (the low frequency component).  

Stock data and FIR filters? Of course. Stock data are
non-starionary. FIR filters compute the weighted avreage
of a finite number of data points. If there is a trend in
the data which is far longer than the FIR filter -- and
such a trend exists in stock data -- this is exactly
what one would expect to see.

> Filtering seems to be happening,
> in that the higher frequencies are being filtered out (the graph looks
> smoother and smoother as the number of days is increased).  The filter
> is acting more like a lowpass filter, rather than bandpass.

Sure. Because the data break with the assumptions of
stationarity etc which the theory of digital filters
is based on.

...
> A few questions:
>     a) how should the sin (k*K)/sin(k) function avoid a divide by zero?
>        (I got around it by starting the array calculation at 1, and
>        calculating the [0] element manually.)

There is no other way.

>     b) is my calculation of the minimum number of coefficients correct?
>        Do I need 2X the number of coeffs in ALL cases, or can I get away
>        with less in low frequencies or high frequencies?

The data don't fit the filter, so do whatever you like.

>     c) is my assumption about the bandwidth size correct (i.e., 8 day
>        period with a 0.8 day bandwidth)?

It is an assumption, i.e. guesswork. Your quess is as good
as anyody else's, until you check whether the assumption
is a good one. Performing such a check may or may not be
possible.

>     d) why doesn't it work? ;-)

If the question still stands by the time you have read to
here, re-read this post.

Rune

On Mar 31, 7:33 pm, info2...@parse.com (Robert Krten) wrote:
>     d) why doesn't it work? ;-)

Maybe your filter does work.  Have you tested it
with data of known spectra?  e.g. random sequences
of various types, pure sine waves, and combinations
thereof?  There are many possible reasons for not
finding an expected output from a bandpass filter;
one is that the filter was broken, another is that
there was no such frequency content in the input.

Ron N. <rhnlogic@yahoo.com> wrote:
> On Mar 31, 7:33 pm, info2...@parse.com (Robert Krten) wrote:
> >     d) why doesn't it work? ;-)

> Maybe your filter does work.  Have you tested it
> with data of known spectra?  e.g. random sequences
> of various types, pure sine waves, and combinations
> thereof?  There are many possible reasons for not
> finding an expected output from a bandpass filter;
> one is that the filter was broken, another is that
> there was no such frequency content in the input.

I have not yet tried it with such testcases -- probably a
very good idea :-)

I *did* manage to fix it, though.  The problem was that
the lowpass filter coeffs I was generating were too high
frequency.  So, where I have:

	f [k] = sin (k * K / N) / sin (k / N);

the value of K that I was using was too large -- this generated
quite a few sine waves as coeffs, and the effect that I need to
generate is more or less a half-wave or less.  Fixed, and it
seems to be operating much better.

I also fixed the shift cos() multiplication to be a sin()
multiplication and that has helped with the output.

I will probably take this function and put it into a library.
At that point, I'll validate it with known data... :-) 

Thanks,
-RK

-- 
Robert Krten, Antique computer collector looking for PDP-8 and PDP-8/S
minicomputers; check out their "good home" at www.parse.com/~museum

Ron N. <rhnlogic@yahoo.com> wrote:
> On Mar 31, 7:33 pm, info2...@parse.com (Robert Krten) wrote:
> >     d) why doesn't it work? ;-)

> Maybe your filter does work.  Have you tested it
> with data of known spectra?  e.g. random sequences
> of various types, pure sine waves, and combinations
> thereof?  There are many possible reasons for not
> finding an expected output from a bandpass filter;
> one is that the filter was broken, another is that
> there was no such frequency content in the input.

I have not yet tried it with such testcases -- probably a
very good idea :-)

I *did* manage to fix it, though.  The problem was that
the lowpass filter coeffs I was generating were too high
frequency.  So, where I have:

	f [k] = sin (k * K / N) / sin (k / N);

the value of K that I was using was too large -- this generated
quite a few sine waves as coeffs, and the effect that I need to
generate is more or less a half-wave or less.  Fixed, and it
seems to be operating much better.

I also fixed the shift cos() multiplication to be a sin()
multiplication and that has helped with the output.

I will probably take this function and put it into a library.
At that point, I'll validate it with known data... :-) 

Thanks,
-RK

-- 
Robert Krten, Antique computer collector looking for PDP-8 and PDP-8/S
minicomputers; check out their "good home" at www.parse.com/~museum

Rune Allnor <allnor@tele.ntnu.no> wrote:
> On 1 Apr, 04:33, info2...@parse.com (Robert Krten) wrote:
> > Hi folks,
> >
> > I'm trying to implement an FIR bandpass filter in C; my reference material
> > is the excellent "Understanding Digital Signal Processing (2nd ed)" by
> > Richard G Lyons.
> >
> > I'm starting with the lowpass filter, and then modifying the parameters
> > for bandpass (described below).
> >
> > I'm calculating the lowpass filter coefficients based on the formula
> > on page 169:
> >
> >         1    sin (PI * k * K / N)
> > h (k) = -- x --------------------
> >         N      sin (PI * k / N)
> >
> > I'm filtering stock data, so my sampling frequency is one sample per day,
> > and I'm interested in bandpass filtering 8 day periods, 16 day, 32 day,
> > and higher, right up to multi-year trends.

> One cycle / year = 1/ 300 or whatever the number of trading
> days per year is... narrow bands. You might need to get
> into advanced techniques to get this one to work.

Fixed.  My "K" value was too high (I had computed it incorrectly).  With
a much lower K, (say 1, or 2) the bandpass that's being generated is
now in line with expectations.  I had been generating a few cycles of
sine wave as coeffs, resulting in a bandpass that was approaching fs/2 :-(

> > The problem I'm running into is that the output of the filter follows the
> > input shape (the low frequency component).  

> Stock data and FIR filters? Of course. Stock data are
> non-starionary. FIR filters compute the weighted avreage
> of a finite number of data points. If there is a trend in
> the data which is far longer than the FIR filter -- and
> such a trend exists in stock data -- this is exactly
> what one would expect to see.

That would apply to the very long sequences, i.e., the "years" samples,
given insufficent length.  But when I have 2 years worth of samples,
and I'm filtering for an 8-day trend, I would expect it to work...

> > Filtering seems to be happening,
> > in that the higher frequencies are being filtered out (the graph looks
> > smoother and smoother as the number of days is increased).  The filter
> > is acting more like a lowpass filter, rather than bandpass.

> Sure. Because the data break with the assumptions of
> stationarity etc which the theory of digital filters
> is based on.

It works now :-)

> > A few questions:
> >     a) how should the sin (k*K)/sin(k) function avoid a divide by zero?
> >        (I got around it by starting the array calculation at 1, and
> >        calculating the [0] element manually.)

> There is no other way.

Ok.  It seemed like a "kludge", but I'm glad that there isn't a more
elegant way.

> >     b) is my calculation of the minimum number of coefficients correct?
> >        Do I need 2X the number of coeffs in ALL cases, or can I get away
> >        with less in low frequencies or high frequencies?

> The data don't fit the filter, so do whatever you like.

Surely some of the data fits the filter model.  Given 2000 samples, and the
fact that I'm looking for 8-sample waves -- this is fine, normal, and
expected, right?

> >     c) is my assumption about the bandwidth size correct (i.e., 8 day
> >        period with a 0.8 day bandwidth)?

> It is an assumption, i.e. guesswork. Your quess is as good
> as anyody else's, until you check whether the assumption
> is a good one. Performing such a check may or may not be
> possible.

The assumption seems ok, given the fix -- which was to 1/x the bandwidth
compution in the code.  Now I'm getting initial lowpass coeffs that are
more like a straight line (i.e., half-wave or less centered at PI/2.)

> >     d) why doesn't it work? ;-)

> If the question still stands by the time you have read to
> here, re-read this post.

Or fix my code :-)

Cheers,
-RK

-- 
Robert Krten, Antique computer collector looking for PDP-8 and PDP-8/S
minicomputers; check out their "good home" at www.parse.com/~museum

Rune Allnor <allnor@tele.ntnu.no> wrote:
> On 1 Apr, 04:33, info2...@parse.com (Robert Krten) wrote:
> > Hi folks,
> >
> > I'm trying to implement an FIR bandpass filter in C; my reference material
> > is the excellent "Understanding Digital Signal Processing (2nd ed)" by
> > Richard G Lyons.
> >
> > I'm starting with the lowpass filter, and then modifying the parameters
> > for bandpass (described below).
> >
> > I'm calculating the lowpass filter coefficients based on the formula
> > on page 169:
> >
> >         1    sin (PI * k * K / N)
> > h (k) = -- x --------------------
> >         N      sin (PI * k / N)
> >
> > I'm filtering stock data, so my sampling frequency is one sample per day,
> > and I'm interested in bandpass filtering 8 day periods, 16 day, 32 day,
> > and higher, right up to multi-year trends.

> One cycle / year = 1/ 300 or whatever the number of trading
> days per year is... narrow bands. You might need to get
> into advanced techniques to get this one to work.

Fixed.  My "K" value was too high (I had computed it incorrectly).  With
a much lower K, (say 1, or 2) the bandpass that's being generated is
now in line with expectations.  I had been generating a few cycles of
sine wave as coeffs, resulting in a bandpass that was approaching fs/2 :-(

> > The problem I'm running into is that the output of the filter follows the
> > input shape (the low frequency component).  

> Stock data and FIR filters? Of course. Stock data are
> non-starionary. FIR filters compute the weighted avreage
> of a finite number of data points. If there is a trend in
> the data which is far longer than the FIR filter -- and
> such a trend exists in stock data -- this is exactly
> what one would expect to see.

That would apply to the very long sequences, i.e., the "years" samples,
given insufficent length.  But when I have 2 years worth of samples,
and I'm filtering for an 8-day trend, I would expect it to work...

> > Filtering seems to be happening,
> > in that the higher frequencies are being filtered out (the graph looks
> > smoother and smoother as the number of days is increased).  The filter
> > is acting more like a lowpass filter, rather than bandpass.

> Sure. Because the data break with the assumptions of
> stationarity etc which the theory of digital filters
> is based on.

It works now :-)

> > A few questions:
> >     a) how should the sin (k*K)/sin(k) function avoid a divide by zero?
> >        (I got around it by starting the array calculation at 1, and
> >        calculating the [0] element manually.)

> There is no other way.

Ok.  It seemed like a "kludge", but I'm glad that there isn't a more
elegant way.

> >     b) is my calculation of the minimum number of coefficients correct?
> >        Do I need 2X the number of coeffs in ALL cases, or can I get away
> >        with less in low frequencies or high frequencies?

> The data don't fit the filter, so do whatever you like.

Surely some of the data fits the filter model.  Given 2000 samples, and the
fact that I'm looking for 8-sample waves -- this is fine, normal, and
expected, right?

> >     c) is my assumption about the bandwidth size correct (i.e., 8 day
> >        period with a 0.8 day bandwidth)?

> It is an assumption, i.e. guesswork. Your quess is as good
> as anyody else's, until you check whether the assumption
> is a good one. Performing such a check may or may not be
> possible.

The assumption seems ok, given the fix -- which was to 1/x the bandwidth
compution in the code.  Now I'm getting initial lowpass coeffs that are
more like a straight line (i.e., half-wave or less centered at PI/2.)

> >     d) why doesn't it work? ;-)

> If the question still stands by the time you have read to
> here, re-read this post.

Or fix my code :-)

Cheers,
-RK

-- 
Robert Krten, Antique computer collector looking for PDP-8 and PDP-8/S
minicomputers; check out their "good home" at www.parse.com/~museum

Rune Allnor <allnor@tele.ntnu.no> wrote:
> On 1 Apr, 04:33, info2...@parse.com (Robert Krten) wrote:
> > Hi folks,
> >
> > I'm trying to implement an FIR bandpass filter in C; my reference material
> > is the excellent "Understanding Digital Signal Processing (2nd ed)" by
> > Richard G Lyons.
> >
> > I'm starting with the lowpass filter, and then modifying the parameters
> > for bandpass (described below).
> >
> > I'm calculating the lowpass filter coefficients based on the formula
> > on page 169:
> >
> >         1    sin (PI * k * K / N)
> > h (k) = -- x --------------------
> >         N      sin (PI * k / N)
> >
> > I'm filtering stock data, so my sampling frequency is one sample per day,
> > and I'm interested in bandpass filtering 8 day periods, 16 day, 32 day,
> > and higher, right up to multi-year trends.

> One cycle / year = 1/ 300 or whatever the number of trading
> days per year is... narrow bands. You might need to get
> into advanced techniques to get this one to work.

Fixed.  My "K" value was too high (I had computed it incorrectly).  With
a much lower K, (say 1, or 2) the bandpass that's being generated is
now in line with expectations.  I had been generating a few cycles of
sine wave as coeffs, resulting in a bandpass that was approaching fs/2 :-(

> > The problem I'm running into is that the output of the filter follows the
> > input shape (the low frequency component).  

> Stock data and FIR filters? Of course. Stock data are
> non-starionary. FIR filters compute the weighted avreage
> of a finite number of data points. If there is a trend in
> the data which is far longer than the FIR filter -- and
> such a trend exists in stock data -- this is exactly
> what one would expect to see.

That would apply to the very long sequences, i.e., the "years" samples,
given insufficent length.  But when I have 2 years worth of samples,
and I'm filtering for an 8-day trend, I would expect it to work...

> > Filtering seems to be happening,
> > in that the higher frequencies are being filtered out (the graph looks
> > smoother and smoother as the number of days is increased).  The filter
> > is acting more like a lowpass filter, rather than bandpass.

> Sure. Because the data break with the assumptions of
> stationarity etc which the theory of digital filters
> is based on.

It works now :-)

> > A few questions:
> >     a) how should the sin (k*K)/sin(k) function avoid a divide by zero?
> >        (I got around it by starting the array calculation at 1, and
> >        calculating the [0] element manually.)

> There is no other way.

Ok.  It seemed like a "kludge", but I'm glad that there isn't a more
elegant way.

> >     b) is my calculation of the minimum number of coefficients correct?
> >        Do I need 2X the number of coeffs in ALL cases, or can I get away
> >        with less in low frequencies or high frequencies?

> The data don't fit the filter, so do whatever you like.

Surely some of the data fits the filter model.  Given 2000 samples, and the
fact that I'm looking for 8-sample waves -- this is fine, normal, and
expected, right?

> >     c) is my assumption about the bandwidth size correct (i.e., 8 day
> >        period with a 0.8 day bandwidth)?

> It is an assumption, i.e. guesswork. Your quess is as good
> as anyody else's, until you check whether the assumption
> is a good one. Performing such a check may or may not be
> possible.

The assumption seems ok, given the fix -- which was to 1/x the bandwidth
compution in the code.  Now I'm getting initial lowpass coeffs that are
more like a straight line (i.e., half-wave or less centered at PI/2.)

> >     d) why doesn't it work? ;-)

> If the question still stands by the time you have read to
> here, re-read this post.

Or fix my code :-)

Cheers,
-RK

-- 
Robert Krten, Antique computer collector looking for PDP-8 and PDP-8/S
minicomputers; check out their "good home" at www.parse.com/~museum

On 2 Apr, 03:20, info2...@parse.com (Robert Krten) wrote:
> Rune Allnor <all...@tele.ntnu.no> wrote:
> > On 1 Apr, 04:33, info2...@parse.com (Robert Krten) wrote:

> > > The problem I'm running into is that the output of the filter follows the
> > > input shape (the low frequency component).  
> > Stock data and FIR filters? Of course. Stock data are
> > non-starionary. FIR filters compute the weighted avreage
> > of a finite number of data points. If there is a trend in
> > the data which is far longer than the FIR filter -- and
> > such a trend exists in stock data -- this is exactly
> > what one would expect to see.
>
> That would apply to the very long sequences, i.e., the "years" samples,
> given insufficent length.  But when I have 2 years worth of samples,
> and I'm filtering for an 8-day trend, I would expect it to work...

You don't get my point: The stock market data don't comply with
the assumtions about data the FIR filter techniques is based on.

The types of signals which are analyzed in the realms of DSP
are assumed to be of type "finite energy" or "finite power".
Such signals are a bit restricted with respect to arbitrary
data sequences in that they have a finite mean value. This
mean value is often - but not always - zero.

Stock market data do not exhibit this limitation. While only
finite amounts of such data have been obsreved, the general
trend is that such data grow with time. So if you want to
analyze these sorts of data, you need to look elsewhere than
DSP for analysis techniques.

Since the FIR filter is a numerical routine, it will produce a
numerical result, and in that sense, it might be said to "work".
However, _a_ numerical result does not mean a _useful_ or
_relevant_ numerical result.

Rune

Rune Allnor <allnor@tele.ntnu.no> wrote:
> On 2 Apr, 03:20, info2...@parse.com (Robert Krten) wrote:
> > Rune Allnor <all...@tele.ntnu.no> wrote:
> > > On 1 Apr, 04:33, info2...@parse.com (Robert Krten) wrote:

> > > > The problem I'm running into is that the output of the filter follows the
> > > > input shape (the low frequency component).  
> > > Stock data and FIR filters? Of course. Stock data are
> > > non-starionary. FIR filters compute the weighted avreage
> > > of a finite number of data points. If there is a trend in
> > > the data which is far longer than the FIR filter -- and
> > > such a trend exists in stock data -- this is exactly
> > > what one would expect to see.
> >
> > That would apply to the very long sequences, i.e., the "years" samples,
> > given insufficent length.  But when I have 2 years worth of samples,
> > and I'm filtering for an 8-day trend, I would expect it to work...

> You don't get my point: The stock market data don't comply with
> the assumtions about data the FIR filter techniques is based on.

> The types of signals which are analyzed in the realms of DSP
> are assumed to be of type "finite energy" or "finite power".
> Such signals are a bit restricted with respect to arbitrary
> data sequences in that they have a finite mean value. This
> mean value is often - but not always - zero.

> Stock market data do not exhibit this limitation. While only
> finite amounts of such data have been obsreved, the general
> trend is that such data grow with time. So if you want to
> analyze these sorts of data, you need to look elsewhere than
> DSP for analysis techniques.

> Since the FIR filter is a numerical routine, it will produce a
> numerical result, and in that sense, it might be said to "work".
> However, _a_ numerical result does not mean a _useful_ or
> _relevant_ numerical result.

Ok.  I get your point now :-)  I will need to do some more research,
obviously :-)

Thanks again,
-RK

-- 
Robert Krten, Antique computer collector looking for PDP-8 and PDP-8/S
minicomputers; check out their "good home" at www.parse.com/~museum