LeafLabs Garden » Topic: Digital filtering

LeafLabs Garden » Topic: Digital filtering http://forums.leaflabs.com/topic.php?id=10396 A place to share, learn, and grow... en-US Fri, 22 Jan 2016 00:22:32 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://forums.leaflabs.com/search.php mlundinse on "Digital filtering" http://forums.leaflabs.com/topic.php?id=10396#post-22972 Mon, 18 Mar 2013 14:11:29 +0000 mlundinse 22972@http://forums.leaflabs.com/ Check the pow call, perhaps it should be pow( resonance + 24/16, 0.5 ) to calculate a square root. Check the values of r and c with a debug println so the filter is stable With array values between 0 and 1024 this corresponds to a very small amplitude, the Maple pwm resolution is 16 bits so scale it up by a factor of 32. As long as the filter is stable the calculated values of the filtered signal, v1, are in the same range as the inputs, so the filtered values can be directly fed to the pwm. int v0,v1; int ci, omrc; float c = 2*3.14*filterval/samplerate; //Small angles, tanx = x, samplerate is 28014hz float r = pow(resonance+24/16, 0.5); //Resonance value lies in the range of 1 - 32ish /* convert filter coefficients to 16.16 fixed points /* Both int and long are 32 bits on arm architecture */ ci = (int)(c*65536); omrc = 65536-(int)(r*c*65536); /* In sample rate loop, Update the phase value, IncrementV2, then read wave table and update filter */ v0 = (omrc*v0 - ci*v1 + ci*32*arrayu[IncrementV2]) >> 16; v1 = (omrc*v1 + ci*v0) >> 16; /* Write output */ pwmWrite(BOARD_ANALOG_OUT, v1); Good luck, I have not tested it or even tried to compile it, but hope it helps. As you have Matlab you can use Matlab to debug the sanity of the algorithms used. Magnus pyrohaz on "Digital filtering" http://forums.leaflabs.com/topic.php?id=10396#post-22971 Mon, 18 Mar 2013 09:53:01 +0000 pyrohaz 22971@http://forums.leaflabs.com/ No problem, thank you for the help so far anyway! The values in arrayu are just the values for a square wave between 1 and 1024: int arrayu[512] = {1024, 1024, 1023, 1024, 1024, 1023, 1024, 1024, 1023, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1023, 1024, 1024, 1023, 1024, 1023, 1024, 1024, 1024, 1024, 1024, 1024, 1023, 1024, 1024, 1024, 1023, 1024, 1024, 1024, 1024, 1024, 1023, 1024, 1024, 1024, 1023, 1024, 1024, 1024, 1024, 1023, 1024, 1023, 1023, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1023, 1023, 1024, 1024, 1024, 1023, 1024, 1024, 1023, 1024, 1024, 1024, 1024, 1023, 1024, 1024, 1024, 1024, 1024, 1023, 1024, 1023, 1024, 1024, 1023, 1023, 1024, 1024, 1024, 1024, 1024, 1023, 1024, 1024, 1023, 1024, 1024, 1023, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1023, 1024, 1024, 1023, 1024, 1024, 1024, 1024, 1023, 1024, 1023, 1023, 1024, 1024, 1024, 1024, 1024, 1023, 1024, 1023, 1023, 1024, 1023, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1023, 1023, 1024, 1024, 1023, 1023, 1024, 1024, 1024, 1023, 1024, 1024, 1023, 1023, 1023, 1023, 1024, 1024, 1024, 1024, 1023, 1024, 1024, 1024, 1024, 1024, 1023, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1023, 1023, 1023, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1023, 1024, 1024, 1024, 1023, 1024, 1024, 1024, 1023, 1024, 1023, 1024, 1023, 1024, 1024, 1024, 1024, 1023, 1023, 1024, 1024, 1023, 1024, 1024, 1023, 1024, 1024, 1023, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1024, 1023, 1023, 1024, 1024, 1024, 1024, 1023, 1024, 1023, 1024, 1024, 1024, 1024, 1024, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; I created them on matlab, hence why a few are 1023 instead but thats not a problem. So sorry to be a pain but would you be able to write me the final piece of code that would be required, all relating v0, v1, r and c? Thanks a bunch! mlundinse on "Digital filtering" http://forums.leaflabs.com/topic.php?id=10396#post-22969 Mon, 18 Mar 2013 02:00:56 +0000 mlundinse 22969@http://forums.leaflabs.com/ yes I was a bit lazy, as you point out, the v1 must also be bitshifted after the multiplications. The values in the arrayu must also be scaled as 16.16 bit fixed point values (is this a sine table ?), and the output v1 scaled to fit the number of bits in the pwm. pyrohaz on "Digital filtering" http://forums.leaflabs.com/topic.php?id=10396#post-22967 Sun, 17 Mar 2013 21:53:31 +0000 pyrohaz 22967@http://forums.leaflabs.com/ gbulmer - unfortunately, I know absolutely nothing about the command line tools! Would they be hard to implement and would they increase the speed of the filter command i'm trying to use? mlundinse - cheers for the help! I tried the code you gave me and it didn't seem to work properly, the output just seemed to produce white noise. Since the code is dependent on both v0 and v1, do I need to do any bitshifts for v1? v0 = (omrc*v0 - (cfp)*v1 + (c)*arrayu[IncrementV2]) >> 16; v1 = (omrc*v1 + (cfp)*v0) >> 16; Would the above work? Sorry for my lack of ingenuity, this is the final step of my project and it would be great if it worked haha. I just realised that as you say scaling by 2^32 because you can rearrange it to (0.1*0.3) * (2^16*2^16) which using simple indices maths turns out as (0.1*0.3*2^32). Cheers so far! mlundinse on "Digital filtering" http://forums.leaflabs.com/topic.php?id=10396#post-22964 Sun, 17 Mar 2013 17:00:53 +0000 mlundinse 22964@http://forums.leaflabs.com/ as long as there is only two factors in the product you can sum them before shifting. r and c are float, dont multiply them with 2^16 until you convert them to fixed point integers. fixed point versions of the filter coefficients: int oneminusrc = 1<<16-(int)(r*c*65536); int cfp = (int)(c*65536) v0 = (oneminusrc*v0 - cfp*v1 + cfp*arrayu[IncrementV1])>>16; Your example, in fixed point arithmetic the real, floating point, value is multiplied by 2^16. When you multily two such numbers the real walue will be scaled by 2^32, but then we shift by 16 bits to get the scaling back to 2^16, as a fixed point integer used in your filter loop. Then afterwards the fixed point value value is divided by floating point 65536.0 to get the real float value back. ((0.1*2^16)*(0.3*2^16))>>16 == 0.03*2^16 gbulmer on "Digital filtering" http://forums.leaflabs.com/topic.php?id=10396#post-22961 Sun, 17 Mar 2013 13:04:28 +0000 gbulmer 22961@http://forums.leaflabs.com/ pyrohaz - are you comfortable with the command-line tools (used under the covers by Maple IDE?) I ask because the <a href="https://launchpad.net/gcc-arm-embedded">ARM Launchpad embedded gcc</a> could either replace the parts in the IDE, or used from the command line. The important point about the Launchpad gcc is it's version 4.7, which supports fixed point, and saturating data types using integer operations <a href="http://gcc.gnu.org/onlinedocs/gcc/Fixed_002dPoint.html" rel="nofollow">http://gcc.gnu.org/onlinedocs/gcc/Fixed_002dPoint.html</a> There is a thread on the forum where someone has used it. I think this is non-trivial, and would take time and effort to get comfortable. However if you're intending to do a lot of coding, it might make the arithmetic expressions easy. pyrohaz on "Digital filtering" http://forums.leaflabs.com/topic.php?id=10396#post-22959 Sun, 17 Mar 2013 11:06:09 +0000 pyrohaz 22959@http://forums.leaflabs.com/ Hey, cheers for the reply! Well the v0 and v1 section is done within an interrupt to write the next sample to the PWM register. Am I right in assuming: filterval = (100+analogRead(filterpot)); c = pow(2, 16)*(6.28*(100+filterval)/SampleRate); r = pow(2, 16)*pow(0.5, ((filterval/128)+30)/16); v0 = ((65535-r*c)*v0 - (c)*v1 + (c)*arrayu[IncrementV1])>>16; v1 = ((65536-r*c)*v1 + (c)*v0)>>16; Or do I need to shift every multiplication by 16? By multiplying the c by 2^16, gives a range from 63114 to 1483 (within range of unsigned int) and multiplying r by 2^16 gives a range of 129818 to 89750 (within the range of a long). Am I correct in assuming that (1-r*c) needs to become (2^16-r*c)? The only reason I really find this hard to grasp is for example 0.1*0.3 in integer math would be ((0.1*2^16)*(0.3*2^16))>>32 Where the floats (0.1 and 0.3) are being multiplied to make them large integers then divided later again to turn them into floats if you get me. Thanks, Harris mlundinse on "Digital filtering" http://forums.leaflabs.com/topic.php?id=10396#post-22958 Sun, 17 Mar 2013 08:20:08 +0000 mlundinse 22958@http://forums.leaflabs.com/ The slow part here is probably the calculation of the filter coefficients. So if you read the filterpot at a lower frequency (10Hz ?) and reclaculate r and c only when the pot values or the resonance have changed, then the main filter loop should run real fast. Use 16.16 or 8.24 fixed point arithmetic. That means the low order 16 or 8 bits of a 32 bit integer are fraction, and the high bits are integer part. - Add numbers as usual - Multiply numbers and then normalize by shifting 16 or 8 bits to the right. - Flot to pixed point, just multiply with the power of 2, 65536 or 256. Regards, Magnus pyrohaz on "Digital filtering" http://forums.leaflabs.com/topic.php?id=10396#post-22955 Sat, 16 Mar 2013 22:08:20 +0000 pyrohaz 22955@http://forums.leaflabs.com/ Hmm after searcing the MusicDSP database, I implemented the 1-RC and C filter. So far it seems to be working well though its being used in a realtime situation and unfortunately seems to be running quite slow. Since its dependent on float math, how can I convert it to use integer math instead? filterval = 1+analogRead(filterpot) //Between 1 and 4096 c = 2*3.14*filterval/samplerate //Small angles, tanx = x, samplerate is 28014hz r = pow(0.5, resonance+24/16)) //Resonance value lies in the range of 1 - 32ish v0 = (1-r*c)*v0 - (c)*v1 + (c)*inputval; v1 = (1-r*c)*v1 + (c)*v0; Output is v1, written to pwm, inputval is a sample in the range of 1 and 1024. Could anyone aid me in converting this to integer math? Thanks, Harris pyrohaz on "Digital filtering" http://forums.leaflabs.com/topic.php?id=10396#post-22951 Sat, 16 Mar 2013 09:10:30 +0000 pyrohaz 22951@http://forums.leaflabs.com/ Hey again guys (sorry for the mass questions!) I'm wondering how I can make a variable low pass filter in code? It needs to be a chebyshev lowpass with a passband ripple of minimum 6dB (to give that resonant peak), cutoff frequency needs to be variable from about 100 up to about 5k or so. Sample rate is 28014Hz and my wavetable consists of values between 0 and 1024. I had a go with the website: <a href="http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html" rel="nofollow">http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html</a> in creating a couple of filters and these work in code but I can't figure out how to vary them! I know its dependent on the output gain, and the two values that preceed yv[0] and yv[1] but I can't seem to find any correlation in how these affect the cutoff frequency. As you can probably tell, i'm not very versed in filter knowledge! Cheers,