Hi All, I'm fairly new to DSP, but have learned a lot since I've started work on it for my senior project at a university. Meaning, please bear with me. For purposes of discussion here, I'm trying to implement IIR filters in the audio frequency range. Accomplishing this has proven to be more difficult than I think it should. I've written Matlab programs to help me along the way. I've written a few programs which calculate the filter coefficients for 2nd and 4th order high pass, bandpass and low pass butterworth filters. These simulate great in matlab. Next step: I've written a matlab program which takes those fractional coefficients, divides by a scalar chosen by the user, and then converts those new fractions to binary and then to hex for actual use on my DSP chip (which is a Freescale DSP56366, by the way. I have the Freescale DSPAUDIOEVM). I simulated with simply dividing all of my coefficients by the largest needed integer to get all coefficients between -1 and 1. This obviously didn't work, as I've learned more about now. I guess to cut to the chase: I'm trying to understand more about fixed point arithmetic. Only tonight after several hours did I get any meaningful results, after reading many posts on these boards. My operation: All coefficients are scaled down by 8 (2^3), and then converted to hex. In my code when I go to store my coefficients, I shift them left by 12 bits. This is effectively Q15 representation (Multiplying by 2^15 and dividing by 2^3 gives 2^12). I then perform my mac operations, and when finished, shift right by 12 bits. Filtering occurs, but not at the correct frequencies (designing for cutoff at 200Hz actually is -3dB at 7.3kHz...), and actually with unexpected gain in some situations. I suppose my longwinded question is, am I on the right path? The last problem I am seeing in my way is this scaling of coefficients, which seems to completely alter filter characteristics in some ways (likely something I am doing wrong). I have read as much as I can find about fixed point arithmetic, but any other resources or opinions would be greatly appreciated. Also, if any of my code is desired just ask. Thanks in advance. - Jim
Fixed Point Arithmetic, Scaling Coefficients
Started by ●February 21, 2008
Reply by ●February 21, 20082008-02-21
Hello Jim, it´s not fairly clear to me what you are doing, but I think you have some fix point arithmetic problem (Your coefficients are too high). So: 1) look in your forward and backward coefficients for biggest coefficient 2) take the log2 of the biggest coefficient and round up to the next bigger integer value = BI 3) divide all coefficients by 2^BI. They should be now between [-1 1[ 4) calculate the fractional representation of the values for your DSP 5) adapt your filter structure with a number of left shifts BI in the forward path of the feedback part of your filter structure to make point 5 more clear (it´s the same for higher order filter): x ___ y>----> b_0---------|+ |------->SHIFT left by BI ------>| |__| | | / \ z^-1 z^-1 / \ | | / \ | |---> b_1----- ---------a_1 <----------- Maybe you have to do the feedback part in double precision (round cycle noise problem). Check that your feedback coefficients have the correct sign - Between the above implementation and the output/input of Matlab functions the "a_" coefficients are toggled in sign. Check that you might have no overflow after the FIR part, or some extension register is capable to manage this in the summation point. Hope that helps, Markus www.two-pi.com
Reply by ●February 21, 20082008-02-21
"jag304" <jag304@lehigh.edu> writes:> [...]Hi Jim, If your frequency response is this far off, you're doing something very wrong somewhere. You should be able to explain exactly what any arithmetic operation does on your DSP using the rules of fixed-point arithmetic. I have documented those rules in a paper, "Fixed Point Arithmetic: An Introduction" http://www.digitalsignallabs.com/fp.pdf You should be able to understand exactly how to convert any coefficient generated by Matlab into a fixed-point representation. You should know what happens in the accumulator when you multiply a coefficient by a data value and sum the result. If taking a look at this documentation or thinking on these things produces more questions, then let us know what they are. -- % Randy Yates % "She's sweet on Wagner-I think she'd die for Beethoven. %% Fuquay-Varina, NC % She love the way Puccini lays down a tune, and %%% 919-577-9882 % Verdi's always creepin' from her room." %%%% <yates@ieee.org> % "Rockaria", *A New World Record*, ELO http://www.digitalsignallabs.com
Reply by ●February 21, 20082008-02-21
Maybe my beautiful ASCII art works with dots instead of blanks: x..................___...............................y>----> b_0---------|+ |------->SHIFT left by BI ------>.|................|__|..........................|.... .|.............. /...\........................z^-1... .z^-1.........../.....\.........................|.... .|............./.......\........................|.... .|---> b_1-----.........---------a_1 <-----------.... For the calculation of your fix point arithmetic values, please use Randys link. Fine Regards, Markus www.two-pi.com
Reply by ●February 21, 20082008-02-21
On Feb 21, 8:25 am, "jag304" <jag...@lehigh.edu> wrote:> Hi All, > > I'm fairly new to DSP, but have learned a lot since I've started work on > it for my senior project at a university. Meaning, please bear with me. > > For purposes of discussion here, I'm trying to implement IIR filters in > the audio frequency range. Accomplishing this has proven to be more > difficult than I think it should. I've written Matlab programs to help me > along the way. I've written a few programs which calculate the filter > coefficients for 2nd and 4th order high pass, bandpass and low pass > butterworth filters. These simulate great in matlab. > > Next step: I've written a matlab program which takes those fractional > coefficients, divides by a scalar chosen by the user, and then converts > those new fractions to binary and then to hex for actual use on my DSP > chip (which is a Freescale DSP56366, by the way. I have the Freescale > DSPAUDIOEVM). > > I simulated with simply dividing all of my coefficients by the largest > needed integer to get all coefficients between -1 and 1. This obviously > didn't work, as I've learned more about now. > > I guess to cut to the chase: I'm trying to understand more about fixed > point arithmetic. Only tonight after several hours did I get any > meaningful results, after reading many posts on these boards. > > My operation: All coefficients are scaled down by 8 (2^3), and then > converted to hex. In my code when I go to store my coefficients, I shift > them left by 12 bits. This is effectively Q15 representation (Multiplying > by 2^15 and dividing by 2^3 gives 2^12). > > I then perform my mac operations, and when finished, shift right by 12 > bits. Filtering occurs, but not at the correct frequencies (designing for > cutoff at 200Hz actually is -3dB at 7.3kHz...), and actually with > unexpected gain in some situations. > > I suppose my longwinded question is, am I on the right path? The last > problem I am seeing in my way is this scaling of coefficients, which seems > to completely alter filter characteristics in some ways (likely something I > am doing wrong). I have read as much as I can find about fixed point > arithmetic, but any other resources or opinions would be greatly > appreciated. > > Also, if any of my code is desired just ask. > Thanks in advance. > > - JimTry using FDATool in Matlab to design filters. You have an option to scale the coefficients. If you use Butterworth you can simplify many of them (they become 1 or -1). If I were you I would scale them between -2 and 2. Cheers
Reply by ●February 21, 20082008-02-21
> >Maybe my beautiful ASCII art works with dots instead of blanks: > >x..................___...............................y >>----> b_0---------|+ |------->SHIFT left by BI ------> >.|................|__|..........................|.... >.|.............. /...\........................z^-1... >.z^-1.........../.....\.........................|.... >.|............./.......\........................|.... >.|---> b_1-----.........---------a_1 <-----------.... > > >For the calculation of your fix point arithmetic values, please useRandys>link. > >Fine Regards, >Markus > >www.two-pi.com >Ok, last try: If you copy and paste the ASCII picture in an text editor (with equally spaced fonts) you will see the correct schematic. Fine regards, Markus www.two-pi.com
Reply by ●February 21, 20082008-02-21
it's been a while since i coded for any Mot 56K or 563xx, but i remember a few things about it. On Feb 21, 6:25 am, "jag304" <jag...@lehigh.edu> wrote:> > My operation: All coefficients are scaled down by 8 (2^3),so you're saying that the natural range of your coefficients are meant to be between -8 and just under +8, right? and the 563xx core in your chip has a natural scaling of -1 to just under +1, right? the Mot assembler knows how to deal with constants (immediate addressing mode or "dc" in memory) that are expressed as fractional, base 10 numbers: -1.0000000 <= the_constant <= 0.99999998> and then > converted to hex. In my code when I go to store my coefficients, I shift > them left by 12 bits.*what* are you shifting left 12 bits? an integer? a fixed point fractional number? for a simple biquad filter, your feedback coefficients will be less than 2 in magnitude, your feedforward coefs could be greater (like 8, which seems like a logical range), but to implement that greater range (of +/- 8), you must shift the output value (in accumulator A or B) left 3 bits before storing the (possibly saturated) accumulator into memory or wherever it goes.> This is effectively Q15 representationbut it's a 24-bit DSP. ???> (Multiplying > by 2^15 and dividing by 2^3 gives 2^12). > > I then perform my mac operations, and when finished, shift right by 12 > bits.uh, i don't think you should be doing that. maybe you should shift *left* by 3 bits. do you not get your output from the most significant word (not the 8 guard bits) in the accumulator (a1 or b1)?> Filtering occurs, but not at the correct frequencies (designing for > cutoff at 200Hz actually is -3dB at 7.3kHz...), and actually with > unexpected gain in some situations. > > I suppose my longwinded question is, am I on the right path? The last > problem I am seeing in my way is this scaling of coefficients, which seems > to completely alter filter characteristics in some ways (likely something I > am doing wrong). I have read as much as I can find about fixed point > arithmetic, but any other resources or opinions would be greatly > appreciated. > > Also, if any of my code is desired just ask.yeah, you better put some of that up. r b-j
Reply by ●February 21, 20082008-02-21
Alrighty then. Here is my MatLab code for converting Fractional Numbers to
Binary and then to Hex. The DSP56366 uses fixed point representation from
2^(0) to 2^(-23) where 2^(0) is the sign bit.
function coef2hex(coef,scaling)
format long
a = zeros(1,24);
if(coef<0) %If negative, first bit is 1
a(1) = 1;
newcoef=-coef/scaling; %Make coefficient positive for processing
else
a(1) = 0; %If positive, first bit is 0
newcoef = coef/scaling;
end
for(i=2:1:24) %Generating binary representation
if(newcoef>=2^(-i+1))
a(i) = 1;
newcoef = newcoef - 2^(-i+1);
end
end
if(a(1)==1)
hitFirst1 = 0;
for(i=24:-1:2) %Starting at the LSB, scan left
if(hitFirst1==1) %If already found first bit,
if(a(i)==1) %invert all following bits to the left
a(i)=0;
else
a(i)=1;
end
end
if(hitFirst1==0 && a(i)==1)
hitFirst1=1;
end
end
end
%Done making binary representation
S = int2str(a);
%Time to go to hex
for(i=4:4:24)
h(i/4) = a(i-3)*8 + a(i-2)*4 + a(i-1)*2 + a(i)*1;
if(h(i/4)<10)
T(i/4) = int2str(h(i/4));
elseif(h(i/4)==10)
T(i/4) = 'A';
elseif(h(i/4)==11)
T(i/4) = 'B';
elseif(h(i/4)==12)
T(i/4) = 'C';
elseif(h(i/4)==13)
T(i/4) = 'D';
elseif(h(i/4)==14)
T(i/4) = 'E';
elseif(h(i/4)==15)
T(i/4) = 'F';
end
end
T
Reply by ●February 21, 20082008-02-21
Next up, one of my programs which computes the filter coefficients based on
cutoff frequency and sampling frequency. I got the coefficient formulae by
using Maple to do a bilinear transformation on a 2nd order butterworth low
pass filter. By defining the sinusoids at the beginning I can test
different filter types and settings.
function secondordlowpass(cutoff,fs,scale)
format long
t = 0:(1/fs):0.025;
x = sin(2*pi*200*t) + sin(2*pi*2000*t);
Wo = 2*pi*cutoff;
a1 = (Wo^2)/(4*fs^2+Wo^2+2*sqrt(2)*Wo*fs)
a2 = 2*a1
a3 = a1
b1 = (-8*fs^2 + 2*Wo^2)/(4*fs^2 + Wo^2 + 2*sqrt(2)*Wo*fs)
b2 = (-2*2^(1/2)*Wo*fs+4*fs^2+Wo^2)/(4*fs^2+Wo^2+2*sqrt(2)*Wo*fs)
%Done finding coefficients
subplot(311),plot(t,x);
title('2kHz Sinusoid Modulated On a 200Hz Sinusoid');
xlabel('Time');
ylabel('Magnitude');
y = zeros(1,length(x));
for(n=3:1:length(x)) %For every sample of the input
y(n) = a1*x(n) + a2*x(n-1) + a3*x(n-2) - b1*y(n-1) - b2*y(n-2);
end
subplot(312),plot(t,y);
title(['Sinusoidal Input Filtered by a 2nd Order Low Pass Butterworth with
a ',num2str(cutoff),...
'Hz Cutoff Frequency']);
xlabel('Time');
ylabel('Magnitude');
coef2hex(a1,scale);
coef2hex(a2,scale);
coef2hex(a3,scale);
coef2hex(-b1,scale);
coef2hex(-b2,scale);
I can keep track of a scaling factor here which is used in converting the
fractional numbers to hex, to be sure there is no overflow.
Reply by ●February 21, 20082008-02-21
On Feb 21, 3:20 pm, "jag304" <jag...@lehigh.edu> wrote:> Next up, one of my programs which computes the filter coefficients based on > cutoff frequency and sampling frequency. I got the coefficient formulae by > using Maple to do a bilinear transformation on a 2nd order butterworth low > pass filter.that's one way to get get 'em.> > coef2hex(a1,scale); > coef2hex(a2,scale); > coef2hex(a3,scale); > coef2hex(-b1,scale); > coef2hex(-b2,scale);why are you doing this at all? don't you know that you can define constants in the 563xx asm language as: org y:$100 a0 dc 0.98765432 a1 dc 0.12345678 ... ??? there is no need anywhere i am aware of to convert your coefficients to integer values, let alone hexadecimal. if you're loading these coefs into the DSP via some CPU that sees everything as integers, then you simply scale the coefficients up by 2^23 (if the range is -1 to +1) or 2^20 (if the range is -8 to +8), round to the nearest integer (MATLAB has the round() function, but if it's in C or someplace without round(), you can add 0.5 and round down with the floor() function). something like this: a0_fixed = round(8388608.0*a0); or a0_fixed = floor(8388608.0*a0 + 0.5); if -1 <= a0 < +1 . (8388608 happens to be 2^23.) as best as i can tell, you're making this so much harder than it need be nor should be. r b-j
