🚀 go-pugleaf

Gerry FitzPatrick (fitzpa@eeel.nist.gov) 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: fuer@siemens.co.at	Gudrunstr. 11
Phone: +43-1-60171-5716		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);
	}
}