NetNews Usenet Archive 1993 #3/NN_1993_3.iso/spool/sci/electron/23208

Text File | 1993-01-22 | 9.8 KB | 312 lines

Newsgroups: sci.electronics
Path: sparky!uunet!news.univie.ac.at!hp4at!siemens.co.at!fuer
From: [e-mail] (Gerhard Fuernkranz)
Subject: Re: Digital Filtering
Sender: [e-mail] (Newssoftware)
Message-ID: <[e-mail]>
Date: Fri, 22 Jan 1993 15:27:02 GMT
References: <fitzpa.727655646@maxwell>
Nntp-Posting-Host: nessie.gud.siemens-austria
Organization: SIEMENS Austria Corp., A-1100 Wien, Gudrunstr. 11
X-Newsreader: Tin 1.1 PL5
Lines: 298

Gerry FitzPatrick ([e-mail]) wrote:
: I need to design a digital filter (either a notch filter or
: a low-pass filter).  I will be analyzing data taken from a
: lecroy 7200 scope (in either ascii or binary).  I would like
: I would like to look at the 60  Hz component of the waveform.
:
: Any ideas on how to write a subroutine that would perform this
: digital filtering on stored data?  Any good books on designing
: digital filters?

Once I have written the following test program
to simulate the behaviour of a digital filter.
It must be linked with -lm.

I hope, this will help you.

--
Gerhard Fuernkranz        SIEMENS Austria Corp.
Email: [e-mail]    Gudrunstr. 11
Phone: [phone]        A-1100 Wien

---------- output produced by this program ----------

The example shows the pulse response of
a 10th order digital filter:

input signal:
Samples:
    0: 0.000000                      *
    1: 0.000000                      *
    2: 0.000000                      *
    3: 0.000000                      *
    4: 1.000000                      |                            *
    5: 1.000000                      |                            *
    6: 1.000000                      |                            *
    7: 1.000000                      |                            *
    8: 1.000000                      |                            *
    9: 1.000000                      |                            *
   10: 1.000000                      |                            *
   11: 1.000000                      |                            *
   12: 1.000000                      |                            *
   13: 1.000000                      |                            *
   14: 1.000000                      |                            *
   15: 1.000000                      |                            *
   16: 1.000000                      |                            *
   17: 1.000000                      |                            *
   18: 1.000000                      |                            *
   19: 1.000000                      |                            *
   20: 1.000000                      |                            *
   21: 1.000000                      |                            *
   22: 1.000000                      |                            *
   23: 1.000000                      |                            *
   24: 1.000000                      |                            *
   25: 1.000000                      |                            *
   26: 1.000000                      |                            *
   27: 1.000000                      |                            *
   28: 1.000000                      |                            *
   29: 1.000000                      |                            *
   30: 1.000000                      |                            *
   31: 1.000000                      |                            *
output signal:
Samples:
    0: 0.000000                      *
    1: 0.000000                      *
    2: 0.000000                      *
    3: 0.000000                      *
    4: 0.000331                      *
    5: 0.004683                      *
    6: 0.030922                      *
    7: 0.126287                      | *
    8: 0.355361                      |       *
    9: 0.725740                      |                *
   10: 1.103899                      |                         *
   11: 1.270497                      |                            *
   12: 1.146250                      |                          *
   13: 0.929901                      |                    *
   14: 0.893190                      |                    *
   15: 1.043636                      |                       *
   16: 1.137073                      |                         *
   17: 1.047580                      |                       *
   18: 0.937792                      |                     *
   19: 0.973822                      |                     *
   20: 1.063945                      |                        *
   21: 1.041176                      |                       *
   22: 0.950166                      |                     *
   23: 0.948518                      |                     *
   24: 1.026586                      |                       *
   25: 1.043214                      |                       *
   26: 0.982888                      |                      *
   27: 0.963995                      |                     *
   28: 1.012875                      |                      *
   29: 1.030241                      |                       *
   30: 0.987176                      |                      *
   31: 0.970688                      |                     *

---------- cut here ----------

#include <math.h>

/*
 * State information for 2nd order digital filter.
 * L1 and L2 a re the currently stored values,
 * C0..C2, D0..D2 are the (digital) filter coefficients.
 */

struct filter {
    double L1, L2;
    double D0, D1, D2;
    double C0, C1, C2;
};

void setup_lp_filter (struct filter *f,
    double OmegaA, double A0, double a1, double b1);
void reset_filter (struct filter *f);
void do_filter (struct filter *f, double *in, double *out, int n);

/*
 * Setup a filter structure for a digital low pass filter
 * of 2nd order.
 * OmagaA is the normalized takeover frequency of the filter.
 * OmegaA = sampling frequency (of the input samples) /
 *        desired takeover frequency of the filter.
 * A0 is the DC gain of the filter (normally 1.0).
 * a1 and b1 are the coefficients of a corresponding
 * 2nd order analog low-pass filter with complex
 * transfer function
 *
 *        A0
 *    ---------------------
 *    1 + a1*P + b1*P^2
 *
 */

void
setup_lp_filter (struct filter *f,
    double OmegaA, double A0, double a1, double b1)
{
    /*
     * Setup anlog filter coefficients.
     * d1, d2, c0 are 0 for a low pass filter.
     * Must be changed for other filter types,
     * e.g. high-pass.
     * The analog (complex) transfer function would be:
     *
     *    d0 + d1*P + d2*P^2
     *    ------------------
     *    c0 + c1*P + c2*P^2
     *
     * Depending on c0..d2, this can be a low-pass,
     * high-pass, band-pass, notch-filter, ...
     */

    double d0 = A0, d1 = 0,  d2 = 0,
           c0 = 1,  c1 = a1, c2 = b1;
    /*
     * Convert analog filter coefficients to digital ones.
     * This should work for all filter types.
     * (Z transformation).
     */

    double l = 1.0/tan(M_PI/OmegaA);
    f->L1 = 0;
    f->L2 = 0;
    f->D0 = (d0 - d1 * l + d2 * l * l) / (c0 + c1 * l + c2 * l * l);
    f->D1 = 2.0 * (d0 - d2 * l * l)    / (c0 + c1 * l + c2 * l * l);
    f->D2 = (d0 + d1 * l + d2 * l * l) / (c0 + c1 * l + c2 * l * l);
    f->C0 = (c0 - c1 * l + c2 * l * l) / (c0 + c1 * l + c2 * l * l);
    f->C1 = 2.0 * (c0 - c2 * l * l)    / (c0 + c1 * l + c2 * l * l);
    f->C2 = 1.0;
}

/*
 * Reset Filter to a defined initial condition.
 */

void
reset_filter (struct filter *f)
{
    f->L1 = 0;
    f->L2 = 0;
}

/*
 * Run an array of input samples in[] thru the filter f
 * and produce an array of output samples out[].
 * N is the number of sampled values in the filter.
 */

void
do_filter (struct filter *f, double *in, double *out, int n)
{
    double L1 = f->L1, L2 = f->L2;
#define IN (*in)
#define OUT (*out)

    while (n-- > 0) {
        OUT = f->D2 * IN + L2;
        L2 = L1 + f->D1 * IN - f->C1 * OUT;
        L1 = f->D0 * IN - f->C0 * OUT;
        in++;
        out++;
    }
    f->L1 = L1;
    f->L2 = L2;
#undef IN
#undef OUT
}

#define        N    32        /* Number of samples */
#define        fsampl    4000.0        /* sampling frequency */
#define        fg    1000.0        /* 3dB takeover frequency
                       of the filter */

/*
 * Filter coefficients for a 10th order butterwort (or is it
 * a 10th order 0.5dB tschebyscheff - I can't remember).
 */

double ai[] = { 6.3648, 1.3582, 0.4822, 0.1994, 0.0563 };
double bi[] = { 18.3695, 4.3453, 1.9440, 1.2520, 1.0263 };

main()
{
    double in[N];
    double in2[N];
    double in3[N];
    double in4[N];
    double in5[N];
    double out[N];
    struct filter f[5];
    int i;

    /*
     * Setup the input signal.
     * Example: step pulse to analyze the
     * pulse response of the filter.
     */

    in[0] = 0;
    in[1] = 0;
    in[2] = 0;
    in[3] = 0;
    for (i = 4; i < N; i++)
        in[i] = 1.0;

    printf("input signal:\n");
    printsamp(in,N);

    /*
     * Set up 10th order low pass filter consisting
     * of 5 2nd order stages.
     */

    for (i = 0; i < 5; i++)
        setup_lp_filter(&f[i],fsampl/fg,1.0,ai[i],bi[i]);

    /*
     * Run the signal thru the 5 filter stages ...
     */

    do_filter(&f[0],in ,in2,N);
    do_filter(&f[1],in2,in3,N);
    do_filter(&f[2],in3,in4,N);
    do_filter(&f[3],in4,in5,N);
    do_filter(&f[4],in5,out,N);

    printf("output signal:\n");
    printsamp(out,N);
}

printsamp (samp, n)
double *samp;
unsigned n;
{
    unsigned i;
    char s[62];
    char nn[32];
    double maxi = 0.0;

    printf ("Samples:\n");

    for (i = 0; i < n; i++)
        if (fabs(samp[i]) > maxi)
            maxi = fabs(samp[i]);
    for (i = 0; i < n; i++) {
        int j;
        memset(s,' ',60);
        s[60] = 0;
        s[30] = '|';
        j = samp[i] / maxi * 30 + 30;
        if (j < 0) j = 0;
        if (j > 59) j = 59;
        sprintf(nn,"%f",samp[i]);
        memcpy(s,nn,strlen(nn));
        s[j] = '*';
        printf("%5d: %s\n",i,s);
    }
}