Maybe someone can put this in laymans terms!
bambino
Posts: 789
I have a spec for a sensor measurement that words like this.
Minimum of 16000 samples per second, through a Low Pass 4-pole butterworth filter having a corner frequency of 1650hz. Digital filtering can be substituted for the butterworth filter.
The 16k sps I have no problem with, I'm actually getting around 56k samples per second. But the digital filtering is giving me a headache. does 1650 hertz mean that i can average my data down to 1650 samples per second?
What is the digital equivalant of a 4-pole butterworth?
A little background: I am taking data from an accelerometer and plotting it in a graph. the actual measurement window is only about 18 ms so I am only dealing with less than 2000 raw samples per graph!
Any comments appreciated, thanks!
Minimum of 16000 samples per second, through a Low Pass 4-pole butterworth filter having a corner frequency of 1650hz. Digital filtering can be substituted for the butterworth filter.
The 16k sps I have no problem with, I'm actually getting around 56k samples per second. But the digital filtering is giving me a headache. does 1650 hertz mean that i can average my data down to 1650 samples per second?
What is the digital equivalant of a 4-pole butterworth?
A little background: I am taking data from an accelerometer and plotting it in a graph. the actual measurement window is only about 18 ms so I am only dealing with less than 2000 raw samples per graph!
Any comments appreciated, thanks!
Comments
Leon
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Amateur radio callsign: G1HSM
Suzuki SV1000S motorcycle
I'd have to port it over to the propeller,Is that code in C or asm!
/*************************************************************************** **************************************************************************** * File: C:/Program Files/mds/dsPIC Filter Design/Butterworth.c * Created by QEDesign Ver 6.3.5 at 15:56:56 Jun 10 2009 * C Code Generator - Version 3.0 **************************************************************************** * Code Fragment to implement filter * * The functions defined in 'qed_filt.c' must be compiled and linked in. * This can be accomplished by either #include "qed_filt.c" * or by separately compiling and linking 'qed_filt.c' * *** following is actual code fragment * extern BiquadSections IIR_Butterworth; * * init_biquad_float (&IIR_Butterworth); // initialize filter structure * * IIR_Butterworth.filter ( x, y, n, &IIR_Butterworth); // x is an array of input samples * // y is an array of output samples * // n is number of samples to process * // &IIR_Butterworth is a pointer to the * // filter structure ***************************************************************************** * This is a complete program which can be compiled and run to test the filter. * To change this to a subroutine only, just add in this program or add globally * in "qed_cgen.h" the line with the definition of DEFINE_SUBROUTINE as follows * #define DEFINE_SUBROUTINE ***************************************************************************** ****************************************************************************/ /* qed_cgen.h contains definitions of filter structures and function prototypes */ #include "qed_cgen.h" /* filter functions are in files 'qed_filt.c' */ float Butterworth_num[noparse][[/noparse] 3] = { 3.372192382813e-001F, /* b[noparse][[/noparse] 1, 0] */ 6.744689941406e-001F, /* b[noparse][[/noparse] 1, 1] */ 3.372192382813e-001F}; /* b[noparse][[/noparse] 1, 2] */ float Butterworth_den[noparse][[/noparse] 3] = { 5.000000000000e-001F, /* a[noparse][[/noparse] 1, 0] */ 6.200256347656e-001F, /* a[noparse][[/noparse] 1, 1] */ 2.289428710938e-001F}; /* a[noparse][[/noparse] 1, 2] */ float Butterworth_m1; float Butterworth_m2; float Butterworth_gain = 1.000000000000e+000F; /* initial gain for cascade realization */ /* also applies to parallel realization */ float Butterworth_pars = 1.000000000000e+000F; /* scale value for parallel sections */ BiquadSections IIR_Butterworth = { 1, /* number of sections */ 0, /* realization method */ /* 0 Cascaded second order sections */ /* 1 Parallel second order sections */ 1, /* quantization: 0 off, 1 on */ 1, /* quantization type */ /* 0 Floating point */ /* 1 Fixed point fractional */ 0, /* realization type for cascaded sections only */ /* 0 Fixed point - Transposed Canonic biquad sections */ /* 1 Fixed point - Canonic biquad sections */ /* 2 Floating point - 4 multiplies */ /* 3 Floating point - 5 multiplies */ 0, /* realization type for parallel sections only */ /* 0 Fixed point - Transposed Canonic biquad sections */ /* 1 Floating point - Transposed Canonic biquad sections */ &Butterworth_gain, /* ptr to gain for cascade/parallel realizations */ &Butterworth_pars, /* ptr to scale value for parallel sections */ Butterworth_num, /* ptr to numerator coefficients */ Butterworth_den, /* ptr to denominator coefficients */ Butterworth_m1, /* ptr to delay line 1 */ Butterworth_m2, /* ptr to delay line 2 */ cas_blkfloat_fm1}; /* ptr to filter routine */ /* call the following function first and normally only once */ /* init_biquad_float (&IIR_Butterworth) */ /* where &IIR_Butterworth is a pointer to the BiquadSections */ /* structure defining the filter */ /* call the following function to filter n samples */ /* IIR_Butterworth.filter (pIn, pOut, int n, &IIR_Butterworth); */ /* where pIn is a pointer to an array or buffer of samples to be filtered */ /* pOut is a pointer to the array of filtered samples output by the filter */ /* n is the number of samples to filter */ /* &IIR_Butterworth is a pointer to the structure defining the filter */ #ifndef DEFINE_SUBROUTINE /* The following main program can be used to test the filter. */ /* input is in file 'in' and the filtered samples are in file 'out' */ /* The input and output files are ascii floating point values */ /* e.g 1.0342 with 1 sample per line */ /* The input files can be created in DSPworks and exported as */ /* ascii floating point or any other system capable of creating */ /* ascii files with floating point values. */ /* The filtered output file can be imported into DSPworks as an ascii */ /* floating point file and an FFT can be run to validate */ /* the frequency response. */ #include "qed_filt.c" #include <stdio.h> #include <string.h> #include <stdlib.h> #include <math.h> #define INSZ1 1000 #define OUTSZ1 1000 static float x[noparse][[/noparse]INSZ1], y[noparse][[/noparse]OUTSZ1]; int main(int argc, char *argv[noparse][[/noparse]]) { int i, in_count, file_status, error; FILE *fin; /* input file of samples */ FILE *fout; /* output file of filtered samples */ fprintf (stderr," ***** start filter test *****\n"); fprintf (stderr," this program accepts 0,1 or 2 command line arguments\n"); fprintf (stderr," the first argument is the filename of the input file\n"); fprintf (stderr," the second argument is the filename of the output file\n"); fprintf (stderr," if there are 0 arguments, input and output is respectively\n"); fprintf (stderr," stdin and stdout\n"); fprintf (stderr," if only one argument is specified, then output is stdout\n"); fprintf (stderr," if input is stdin rather than a file, then fscanf expects input\n"); fprintf (stderr," from the console which may be piped in or entered directly\n"); fin = stdin; fout = stdout; error = 0; if (argc == 1) { fprintf(stderr," ***** waiting for input *****\n"); } if (argc >= 2) { fin = fopen (argv, "r"); if (fin == NULL) { fprintf(stderr,"\n error - Cannot open file %s for input\n", argv); error = 1; } } if (argc >= 3) { fout = fopen (argv, "w"); if (fout == NULL) { fprintf(stderr,"\n error - Cannot open file %s for output\n", argv); error = 1; } } if (error) { fprintf(stderr," ***** end filter test *****\n"); return(0); } init_biquad_float (&IIR_Butterworth); do { /* get input samples */ for (in_count = 0; in_count < INSZ1; in_count++) { file_status = fscanf(fin,"%f",&x[noparse][[/noparse]in_count]); if (file_status != 1) break; } /* filter samples */ if (in_count == 0) break; IIR_Butterworth.filter( x, y, in_count, &IIR_Butterworth); for (i = 0; i < in_count; i++) fprintf (fout,"%f\n",y[i]); } while (file_status == 1); fclose (fin); fclose (fout); fprintf(stderr," ***** end filter test *****\n"); return(1); } #endif [/i]Here is the filter spec:
FILTER SPECIFICATION FILE
FILTER TYPE:LOW PASS 1H
PASSBAND RIPPLE IN -dB -.2000E+01
STOPBAND RIPPLE IN -dB -.3000E+01
PASSBAND CUTOFF FREQUENCY 0.400000E+04 HERTZ
STOPBAND CUTOFF FREQUENCY 0.410000E+04 HERTZ
SAMPLING FREQUENCY 0.100000E+05 HERTZ
The software is fairly expensive.
Leon
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Amateur radio callsign: G1HSM
Suzuki SV1000S motorcycle
Leon
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Amateur radio callsign: G1HSM
Suzuki SV1000S motorcycle
thanks
Post Edited (bambino) : 6/10/2009 3:08:54 PM GMT
IIR_Butterworth.filter( x, y, in_count, &IIR_Butterworth);
for (i = 0; i < in_count; i++)
fprintf (fout,"%f\n",y);
I assume the rest is bundled up in the call to the IIR_Butterworth.filter() routine
>·does 1650 hertz mean that i can average my data down to 1650 samples per second?
Bambino,
corner frequency has to do with the filter's roll-off point - if I recall correctly, it's where the tangent of·the flat-response side and·of the slope of the roll-off curve meet - they "corner" at that point of intersection.
> Low Pass 4-pole butterworth filter,· corner frequency of 1650hz.
#define NZEROS 4 #define NPOLES 4 #define GAIN 1.867017256e+02 static float xv[noparse][[/noparse]NZEROS+1], yv[noparse][[/noparse]NPOLES+1]; static void filterloop() { for (;;) { xv[noparse][[/noparse]0] = xv[noparse][[/noparse]1]; xv[noparse][[/noparse]1] = xv[noparse][[/noparse]2]; xv[noparse][[/noparse]2] = xv[noparse][[/noparse]3]; xv[noparse][[/noparse]3] = xv[noparse][[/noparse]4]; xv[noparse][[/noparse]4] = [i]next input value[/i] / GAIN; yv[noparse][[/noparse]0] = yv[noparse][[/noparse]1]; yv[noparse][[/noparse]1] = yv[noparse][[/noparse]2]; yv[noparse][[/noparse]2] = yv[noparse][[/noparse]3]; yv[noparse][[/noparse]3] = yv[noparse][[/noparse]4]; yv[noparse][[/noparse]4] = (xv[noparse][[/noparse]0] + xv[noparse][[/noparse]4]) + 4 * (xv[noparse][[/noparse]1] + xv[noparse][[/noparse]3]) + 6 * xv[noparse][[/noparse]2] + ( -0.1774344076 * yv[noparse][[/noparse]0]) + ( 1.0075092500 * yv[noparse][[/noparse]1]) + ( -2.2349793544 * yv[noparse][[/noparse]2]) + ( 2.3192063218 * yv[noparse][[/noparse]3]); [i]next output value[/i] = yv[noparse][[/noparse]4]; } }I can't vouch for the accuracy as I've not had time to test this. It was generated here:
http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html
(Yeah, the math behind digital filtering and Fourier Transforms in general require many bottles of Excedrin Extra Strength [noparse]:)[/noparse]
HTH
- Howard in Florida
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Got Electrons?
OK, My C is pretty rusty, but stay tuned and I'll try to post some spin code ported from the filterloop routine.
It may be friday before I'll have time to finish it though.
Thanks for the link. It want help me make a full blown IIR object since I don't know how it arrived at the gain value, but it will let me carry on with my app. Thanks again.
One thing though! Does yv[noparse][[/noparse]4] have valid date output before the 4th loop through the filter or should those samples be dropped?
Post Edited (bambino) : 6/11/2009 12:37:39 PM GMT
···················· + ( -0.1774344076 * yv[noparse][[/noparse]0]) + (· 1.0075092500 * yv[noparse][[/noparse]1])
···················· + ( -2.2349793544 * yv[noparse][[/noparse]2]) + (· 2.3192063218 * yv[noparse][[/noparse]3]);
These fixed numbers here? are they decided from the gain?
What I really need to know is will these numbers be different if I ran different numbers through the calculator you posted the link to? I probably won't, but it would help me understand the routine more!
RE yv and valid results - it's very common that you have to let the filter run a few cycles to preload the values, make the transforms, and get good output. Roughly, this initial lag is usually negligible. As you can see from the code, values are recycled. Generally, this is called 'convolution' --- so yes, you do have to ignore some of the very initial results as they are undefined --- I just pick a number slightly larger than the whole convolution cycle, and ignore those results. Just glancing at that code, it might be four or five samples... but you might want to test that.
RE fixed values: try running different parms through that filter generator, and you'll see how it works - it's pretty straight forward. Yes, those fixed values *are* a direct result of the filter generating parameters.
- H
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Got Electrons?
www.leonheller.com/Designs/Butterworth.gif
Scilab (free) can be used for designing both FIR and IIR filters, and will simulate them.
Leon
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Amateur radio callsign: G1HSM
Suzuki SV1000S motorcycle
@CRP That's great, just adjust the fixed values like I would a PID routine, I presume. Sort of trial and error.
Leon
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Amateur radio callsign: G1HSM
Suzuki SV1000S motorcycle
Right, like most things [noparse]:)[/noparse]
Follow Leon's lead - he understands these things better than I --- and you both have the tools at hand.
cheers
- Howard
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Got Electrons?
I think I've got enough to get started, maybe my honey do list is short this evening!
After thinking about it a while, I've decided to put the filter on the PC end of the App. Not Permenantly, but at least until the initial tests get done. I've got room for the floating point routines and extra cogs available, but the program in the propeller is the most bug free at the moment, so I'll leave it along and write the filter in VB.
Thanks again to you two for all your help!