|
|
View previous topic :: View next topic |
Author |
Message |
cbarberis
Joined: 01 Oct 2003 Posts: 172 Location: Punta Gorda, Florida USA
|
Problem with integer math on digital filter |
Posted: Tue Jul 09, 2019 10:20 am |
|
|
As usual I come to the forum when I’m lost and confused. I have been working on digital filtering for the last few days just to try some new thing, not that I have not used digital filters in the past.
Actually, I was using the free version of filter design software from Microchip but I found it did not work correctly after my translations and decided to look for something else. I came across a great free filter design package on the web called WinFilter ( http://www.winfilter.20m.com ) which is very nice and seems to support most of the filter types that I commonly use.
So I put it to the test doing various IIR filter topologies using floats, all of them worked just fine as I expected. This filter design package also allows to use integer filter math for both IIR and FIR and this is where things went down hill, it just did not work even though there is no obvious reason but it seems to be a numerical issue when the data is put through the filter you just get garbage back typically zeros or ones when you are feeding the filter a 12bit quantity.
This program uses signed int16 for input and output and also uses a signed 32 bit for the sum accumulator but somewhere just before it returns the filtered values it breaks. I am using the latest Microchip Explorer 16/32 development board with a dsPIC33EP256GP506 running at 48MHz and unfortunately the CCS ICD-U80 for some reason does not work on debug mode with this board.
The way I have it set-up is to take a signal from the function generator that is biased properly so my ADC will work with such input. I sample this signal which is typically a sinusoid between 10Hz to 1500Hz, my sampling is done through a timer interrupt at a rate of 5KHz, so I have 200 uS to do both my sampling, filtering and DAC (MCP4921) process and my overhead timing constrain does not seem to be a problem at all for the exception that the higher the frequency I sample the less the number of points and eventually I get distortion and aliasing.
To prove that my hardware was working OK I have done a simple ADC to DAC loop to verify that whatever comes in from the ADC goes to the DAC and I am seeing perfect waveforms up to ~5KHz above that the constrain is due to the DAC overhead latency but I am getting excellent data from analog to digital and back to analog. So I will include some of my code and the code for the filter in question. Just to clarify, there is absolutely no problem with the integer filter constant values to prove this I wrote a simple C test program using Visual Studio.
In this program I basically read an array 50 samples of my ADC captured data at a Fs=5000Hz for a 400Hz sinusoid I then feed each and everyone of the samples to the filter routine which I have included on this test program as a subroutine. The result of this test basically proves that there is nothing wrong with the filter coefficients or the filter math as I am seeing good values coming back from the filter.
However, the same does not happen when the filtered data comes back from the same routine on the actual hardware platform, I am seeing incredibly low numbers returning from the filter data on the passband frequency, these return filtered values also do not correlate to the input values they appear to be non-monotonic.
BTW the float routines are working fine so I have not included them here, also for debug purposes I have created a couple of buffers that captures 50 samples of data from the ADC and the same 50 points of data returned from the filter, this only happens when I press a switch which spits it out of the serial port, otherwise this has no effect on data sampling. I am sure the cause of this problem is something stupid that I am just not seeing . Thank you all for advice.
My code:
Code: |
#include <33EP256GP506.h>
#device ICSP=1
#device ADC=12
#FUSES NODEBUG
#FUSES NOWDT //No Watch Dog Timer
#FUSES CKSFSM //Clock Switching is enabled, fail Safe clock monitor is enabled
#use delay(crystal=8000000,clock= 48000000,pll_wait)
#ifdef INTS
signed int adc_samples[50];
signed int dac_values[50];
signed int adclevel;
signed int filteredval;
#endif
#INT_TIMER1
void timer1_isr(void)
{
#ifdef INTS
//DEBUGPULSE;
adclevel = (signed int)read_adc();
#endif
samplecount++;
if(samplecount >50)
samplecount = 0;
SampleFlag = 1;
}
void Init(void) {
TRISA4 =0;
TRISF1 =0;
TRISC3 =0;
TRISB2 =1;
setup_timer1(TMR_INTERNAL | TMR_DIV_BY_1, 4800); //4800 = 200uS interrupt
setup_adc_ports(sAN0 | sAN2, VSS_VDD);
setup_adc(ADC_CLOCK_DIV_4 | ADC_TAD_MUL_4);
set_adc_channel(2);
init_dac();
enable_interrupts(INT_RDA);
enable_interrupts(INTR_GLOBAL);
}
void main()
{
unsigned int i,n=0;
Init();
enable_interrupts(INT_TIMER1); // Sampling interrupt
while(TRUE)
{
// DEBUGPULSE; // just to get an idea of open loop sampling rate measured out of port B1
#ifdef INTS // filter routine for integer filter math
//DEBUGPULSE;
//adclevel = read_adc();
if( SampleFlag ){
//DEBUGPULSE2;
filteredval = iir(adclevel);
write_dac(filteredval);
adc_samples[n] = adclevel;
//write_dac(adclevel); //only used to verify ADC direct to DAC without going through filter
// DEBUGPULSE2;
dac_values[n] = filteredval;
n++;
if(n >= MAXCAPTURE)
n=0;
CalculateFlag = 0; }
if(SW3DOWN) // Wen you press SW3 on the board
//Prints out 50 0f the incoming ADC values and their corresponding values to the DAC
{ delay_ms(15);
printf("ADCvalue Cnts: DACvalue Cnts \n\r");
for(i=0; i< MAXCAPTURE; i++)
printf(" %ld %ld \n\r",adc_samples[i],dac_values[i]);
}
#endif
}
}
|
Filter Code:
Code: | **************************************************************
Filter type: Band Pass
Filter model: Chebyshev
Filter order: 6
Sampling Frequency: 2000 Hz
Fc1 and Fc2 Frequencies: 250.000000 Hz and 950.000000 Hz
Pass band Ripple: 1.000000 dB
Coefficents Quantization: 16-bit
Z domain Zeros
z = -1.000000 + j 0.000000
z = -1.000000 + j 0.000000
z = -1.000000 + j 0.000000
z = -1.000000 + j 0.000000
z = -1.000000 + j 0.000000
z = -1.000000 + j 0.000000
z = 1.000000 + j 0.000000
z = 1.000000 + j 0.000000
z = 1.000000 + j 0.000000
z = 1.000000 + j 0.000000
z = 1.000000 + j 0.000000
z = 1.000000 + j 0.000000
Z domain Poles
z = -0.770187 + j -0.294026
z = -0.770187 + j 0.294026
z = -0.948519 + j -0.160836
z = -0.948519 + j 0.160836
z = -0.058502 + j -0.589541
z = -0.058502 + j 0.589541
z = 0.463750 + j 0.701452
z = 0.463750 + j -0.701452
z = -0.969485 + j -0.183635
z = -0.969485 + j 0.183635
z = 0.676743 + j -0.679985
z = 0.676743 + j 0.679985
***************************************************************/
#define NCoef 12
#define DCgain 4
int16 iir(int16 NewSample) {
int16 ACoef[NCoef+1] = {
1286,
0,
-7716,
0,
19290,
0,
-25721,
0,
19290,
0,
-7716,
0,
1286
};
int16 BCoef[NCoef+1] = {
4096,
13158,
14330,
8972,
14369,
21201,
16141,
10909,
10055,
7533,
4426,
2211,
573
};
static signed int32 y[NCoef+1]; //output samples
//Warning!!!!!! This variable should be signed (input sample width + Coefs width + 12 )-bit width to avoid saturation.
static int16 x[NCoef+1]; //input samples
int n;
//shift the old samples
for(n=NCoef; n>0; n--) {
x[n] = x[n-1];
y[n] = y[n-1];
}
//Calculate the new output
x[0] = NewSample;
y[0] = ACoef[0] * x[0];
for(n=1; n<=NCoef; n++)
y[0] += ACoef[n] * x[n] - BCoef[n] * y[n];
y[0] /= BCoef[0];
return y[0] / DCgain;
}
|
|
|
|
newguy
Joined: 24 Jun 2004 Posts: 1908
|
|
Posted: Tue Jul 09, 2019 10:30 am |
|
|
Just for s & g, force the order of precedence here:
Code: | y[0] += (ACoef[n] * x[n]) - (BCoef[n] * y[n]); |
|
|
|
cbarberis
Joined: 01 Oct 2003 Posts: 172 Location: Punta Gorda, Florida USA
|
|
Posted: Tue Jul 09, 2019 10:58 am |
|
|
Good try !! But it did not have any effect. BTW I mentioned that by using integer math I would get numbers from the filter that appeared saturated or zeros and that was very noticeable with the HPF filter topology but the BPF filter gives me back garbage data it does not appear to be saturated data. Unfortunately I have no idea as to how to post the images of my scope captures so it can be seeing here, I tried using the Img button but I think it is looking for a URL where I would have to post the images. |
|
|
newguy
Joined: 24 Jun 2004 Posts: 1908
|
|
Posted: Tue Jul 09, 2019 11:24 am |
|
|
Your passband gain is about -20dB. If you put a 5Vp-p sine through your filter, your output in the passband would be about 50mVp-p.
Use 32 bit integers in your filtering routine and scale all inputs to the routine by 128 (<< 7). You should get something out of your DAC then *but* the extra time dealt with using int32's might screw up the loop timing. |
|
|
jeremiah
Joined: 20 Jul 2010 Posts: 1349
|
|
Posted: Tue Jul 09, 2019 11:24 am |
|
|
Since multiplies happen before subtractions, you might try making the x[n] array int32 instead of int16 to see if you are potentially having any int16*int16 multiplication overflows? Particularly from the
Code: | y[0] = ACoef[0] * x[0]; |
and
parts
You might also monitor your y[0] value to make sure you aren't underflowing it by subtracting too large a number multiple times in that for loop. I don't know how big your results from the Bcoeff * y[n] calc get before the subtraction. |
|
|
cbarberis
Joined: 01 Oct 2003 Posts: 172 Location: Punta Gorda, Florida USA
|
|
Posted: Tue Jul 09, 2019 12:04 pm |
|
|
As mentioned before I have taken that exact same filter routine and fed it 50 actual samples obtained at the 5KHz sampling data of course this was using Visual C the output to the filter went into a 50 point buffer which I later inserted into my CCS code and run the program using the same sampling interval feeding the filter output buffer I got from the C test program into the DAC and I obtained a perfectly clear sinusoid.
I realize that integer math computations done on a PC differ quite a bit from an embedded controller. I also suspect that the problem is with the numerical handling, not sure if it is the PIC or the CCS compiler but in any event I will try all your suggestions like declaring the x[n] array as a 32 bit signed int.
I am also seeing something interesting like the filter data is oscillating and being modulated by the real signal, when I take the real signal input from the ADC away the filter keeps oscillating. However, If I start the run with no signal input into the filter routine I get a nice flat line which I should but once a signal is fed into the filter routine the oscillation starts and never goes away.
I know for a fact that all the filter coefficients are calculated from poles and zeros which lie within the Z unit circle so there should be no instability unless numerical value is corrupted.
BTW I did try declaring the x[n] array as a 32 bit signed int. did NOT help. |
|
|
newguy
Joined: 24 Jun 2004 Posts: 1908
|
|
Posted: Tue Jul 09, 2019 1:13 pm |
|
|
Take those same 50 samples, feed them into your project serially, and have a look at the output of the filtering routine (again, serially). Compare with the visual C program's output.
Regarding filter oscillation once fed a real signal: it's an IIR filter. It will do that. |
|
|
cbarberis
Joined: 01 Oct 2003 Posts: 172 Location: Punta Gorda, Florida USA
|
|
Posted: Tue Jul 09, 2019 1:25 pm |
|
|
I actually fed the 50 ADC samples that I used on my C test program to the filter and I got totally different numbers out of the filter which did not make sense compared to the ones I got from my test program. To make things more confusing I noticed that the B coefficients shown on the filter code do not have any negative values which seems odd as the denominator should contain some negative values in most cases they alternate.
So I took the B coefficients as shown and inserted a - (negative value) on very other value of course there was no science to this just guessing in the dark, but any way I got the filter to appear to work, the only problem is that when I take the input signal into the filter away the filter oscillates a very nice sinusoid at ~80Hz does this give anyone any hints?
I am contacting the fellow that wrote this program, I doubt that this never worked but it seems weird. |
|
|
|
|
You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|