From news.kth.se!eru.mt.luth.se!bloom-beacon.mit.edu!newsxfer.itd.umich.edu!news.uoregon.edu!vixen.cso.uiuc.edu!uwm.edu!chi-news.cic.net!newsfeed.internetmci.com!news.sprintlink.net!ddi2.digital.net!picard.csihq.com!jday Thu Nov  2 15:32:26 1995
Path: news.kth.se!eru.mt.luth.se!bloom-beacon.mit.edu!newsxfer.itd.umich.edu!news.uoregon.edu!vixen.cso.uiuc.edu!uwm.edu!chi-news.cic.net!newsfeed.internetmci.com!news.sprintlink.net!ddi2.digital.net!picard.csihq.com!jday
From: jday@csihq.com (John Day)
Newsgroups: comp.dsp
Subject: Re: Simple filter algorithms needed
Date: Sun, 29 Oct 95 17:26:37 GMT
Organization: Computer Science Innovations Inc.
Lines: 101
Distribution: inet
Message-ID: <4709lm$8sg@picard.csihq.com>
References: <46lmb4$hjl@erinews.ericsson.se>
NNTP-Posting-Host: johnday.csihq.com
X-Newsreader: News Xpress Version 1.0 Beta #4

>Can you give me a pointer or post a simple LPF algorithm in C where
>the cutoff frequency can be set.

Hi,

I think I can help you and in return please help me by testing this
LPF/BPF/HPF routine which I have written but never got around to
testing or using.

First, some theory.

All you need to know is how to design an LPF because the BPF and 
HPF can be derived from a LPF by a simple relation I will show you
in a moment.

We want our LPF to have an ideal, 'brickwall' response, which can 
be charactized (in the freq domain) as multiplying by one up to
the cutoff and by zero beyond that up to infinity. 

Recall that multiplication in the freq domain is equivalent to
convolution in time domain (Convolution Theorem). So we can
design our filter in the domain by convolving with the sinc
function which is the inverse transform of the brick wall:

   lpf(n) = sin(fc*(n - n0))/ (pi*(n-n0))

where fc is the cutoff and n0 is the median index. lpf(n) is
the impulse response for your LPF and can be used directly
by convolving it with your time series. Unfortunately, (you
knew there was catch) your filter will need an infinite number
of taps to retain its ideal, brick-wall behaviour. It will 
work with a small number of taps (say 16) but you'll get a
lot of 'leakage'. But you can make it a long as you want.
With a hundred taps or so you can actually get very close
to a brick wall response!

Now to get a BPF just shift the passband by cosine multiplication:
    
    bpf(n) = 2*lpf(n)*cos(wc*n) 
    where wc is the wave number of the passband center.

To get a HPF is even simpler:

    hpf(n) = (-1)**n * lp
    i.e. just flip the sign of the odd coefficients!

All of this is explained very clearly by Prof. J. Cadzow in his
book Foundation of Digital Signal Processing and Data Analysis, 
Chap 5, 'Theory of Frequency-Discrimination Filters'

Now here's some code in rough form. Fix it up, get it working (then
send a copy back to me. Thanks!) ;)

/*------------ Construct ideal lowpass,bandpass or hipass filter ---------*/

vect *vfilt(int ftype, int ntaps, double cutoff ,double centerf)
{
    int    i;
    double n=(double)ntaps,
           nc=(int32)(cutoff*n/PI),
           n0=(n-1.0)/2.0;
    double *filt=malloc(ntaps);

    for (i=0; i<ntaps; i++)
    {
#if 0  // one maybe both expressions are incorrect!
        filt[i] = sin((PI*(2.0*nc+1.0)*(i-n0))/n)/
                        (n*sin((PI*(i-n0))/n));
#else
        filt[i] = sin(cutoff*(i-n0))/(PI*(i-n0));
#endif
        switch (ftype)
        {
            case FLT_HIPASS:
                if (ODD(i))
                    filt[i] *= -1;
                 break;
            case FLT_BANDPASS:
                filt[i] *= 2*cos(i*centerf);
                break;
            case FLT_LOPASS:
            default:
                break;
        } /* endswitch */
    } /* endfor */

    return filt;

} /* end vfilt() */


+---------------------------------------------------------------------------+
|  John Day
|  Computer Science Innovations,Inc
|  Principal Engineer PHONE: (407) 676-2923 ext:410
|  Melbourne, Fl        FAX: (407) 676-2355
|  WWW:   http://www.csihq.com
|  EMAIL: jday@csihq.com
|
|  'God created Integers, all else is the work of Man' -Kronecker
+---------------------------------------------------------------------------+