15
$\begingroup$

I've designed a very simple low-pass Butterworth filter using Matlab. The following code snippet demonstrates what I've done.

fs = 2.1e6;
flow = 44 * 1000;
fNorm =  flow / (fs / 2);
[b,a] = butter(10, fNorm, 'low');

In [b,a] are stored the filter coefficients. I would like to obtain [b,a] as integers so that I can use an online HDL code generator to generate code in Verilog.

The Matlab [b,a] values seem to be too small to use with the online code generator (the server-side Perl script refuses to generate code with the coefficients), and I am wondering if it would be possible to obtain [b,a] in a form that can be used as a proper input.

The a coefficients that I get in Matlab are:

1.0000
-9.1585
37.7780
-92.4225
148.5066
-163.7596
125.5009
-66.0030
22.7969
-4.6694
0.4307

The b coefficients that I get in Matlab are:

1.0167e-012
1.0167e-011
4.5752e-011
1.2201e-010
2.1351e-010
2.5621e-010
2.1351e-010
1.2201e-010
4.5752e-011
1.0167e-011
1.0167e-012

Using the online generator, I would like to design a filter with a 12-bit bitwidth and I or II filter form. I don't know what is meant by the "fractional bits" at the above link.

Running the code generator (http://www.spiral.net/hardware/filter.html) with the [b,a] coefficients listed above, with fractional bits set at 20 and a bitwidth of 12, I receive the following run error:

Integer A constants: 1048576 -9603383 39613104 -96912015 155720456 -171714386 131597231 -69209161 23904282 -4896220 451621
Integer B constants: 0 0 0 0 0 0 0 0 0 0 0

Error: constants wider than 26 bits are not allowed, offending constant = -69209161, effective bitwidth = 7 mantissa + 20 fractional = 27 total.

An error has occurred - please revise the input parameters. 

How might I change my design so that this error does not occur?

UPDATE: Using Matlab to generate a 6th-order Butterworth filter, I get the following coefficients:

For a:

1.0000
-5.4914
12.5848
-15.4051
10.6225
-3.9118
0.6010 

for b:

0.0064e-005
0.0382e-005
0.0954e-005
0.1272e-005
0.0954e-005
0.0382e-005
0.0064e-005

Running the online code generator (http://www.spiral.net/hardware/filter.html), I now receive the following error (with fractional bits as 8 and bitwidth of 20):

./iirGen.pl -A 256  '-1405' '3221' '-3943' '2719' '-1001' '153' -B  '0' '0' '0' '0' '0' '0' '0' -moduleName acm_filter -fractionalBits 8 -bitWidth 20 -inData inData  -inReg   -outReg  -outData outData -clk clk -reset reset -reset_edge negedge -filterForm 1  -debug  -outFile ../outputs/filter_1330617505.v 2>&1 
At least 1 non-zero-valued constant is required.  Please check the inputs and try again.

Perhaps the b-coefficients are too small, or perhaps the code generator (http://www.spiral.net/hardware/filter.html) wants the [b,a] in another format?

UPDATE:

Perhaps what I need to do is scale the [b,a] coefficients by the number of fractional bits to obtain the coefficients as integers.

a .* 2^12
b .* 2^12

However, I still think that the b coefficients are extremely small. What am I doing wrong here?

Perhaps another type of filter (or filter design method) would be more suitable? Could anyone make a suggestion?

UPDATE: As suggested by Jason R and Christopher Felton in the comments below, an SOS filter would be more suitable. I've now written some Matlab code to obtain an SOS filter.

fs = 2.1e6;
flow = 44 * 1000;      
fNorm =  flow / (fs / 2);
[A,B,C,D] = butter(10, fNorm, 'low');
[sos,g] = ss2sos(A,B,C,D);

The SOS matrix that I get is:

1.0000    3.4724    3.1253    1.0000   -1.7551    0.7705
1.0000    2.5057    1.9919    1.0000   -1.7751    0.7906
1.0000    1.6873    1.0267    1.0000   -1.8143    0.8301
1.0000    1.2550    0.5137    1.0000   -1.8712    0.8875
1.0000    1.0795    0.3046    1.0000   -1.9428    0.9598

Is it possible to still use the Verilog code generation tool (http://www.spiral.net/hardware/filter.html) to implement this SOS filter, or should I simply write the Verilog by hand? Is a good reference available?

I would wonder if an FIR filter would be better to use in this situation.

MOREOVER: Recursive IIR filters can be implemented using integer math by expressing coefficients as fractions. (See Smith's excellent DSP signal processing book for further details: http://www.dspguide.com/ch19/5.htm)

The following Matlab program converts Butterworth filter coefficients into fractional parts using the Matlab rat() function. Then as mentioned in the comments, second order sections can be used to numerically implement the filter (http://en.wikipedia.org/wiki/Digital_biquad_filter).

% variables
% variables
fs = 2.1e6;                     % sampling frequency           
flow = 44 * 1000;               % lowpass filter


% pre-calculations
fNorm =  flow / (fs / 2);       % normalized freq for lowpass filter

% uncomment this to look at the coefficients in fvtool
% compute [b,a] coefficients
% [b,a] = butter(7, fNorm, 'low');
% fvtool(b,a)  

% compute SOS coefficients (7th order filter)
[z,p,k] = butter(7, fNorm, 'low');

% NOTE that we might have to scale things to make sure
% that everything works out well (see zp2sos help for 'up' and 'inf' options)
sos = zp2sos(z,p,k, 'up', 'inf'); 
[n,d] = rat(sos); 
sos_check = n ./ d;  % this should be the same as SOS matrix

% by here, n is the numerator and d is the denominator coefficients
% as an example, write the the coefficients into a C code header file
% for prototyping the implementation

 % write the numerator and denominator matices into a file
[rownum, colnum] = size(n);  % d should be the same
sections = rownum;           % the number of sections is the same as the number of rows
fid = fopen('IIR_coeff.h', 'w');

fprintf(fid, '#ifndef IIR_COEFF_H\n');
fprintf(fid, '#define IIR_COEFF_H\n\n\n');
for i = 1:rownum
   for j = 1:colnum

       if(j <= 3)  % b coefficients
            bn = ['b' num2str(j-1) num2str(i) 'n' ' = ' num2str(n(i,j))];
            bd = ['b' num2str(j-1) num2str(i) 'd' ' = ' num2str(d(i,j))];
            fprintf(fid, 'const int32_t %s;\n', bn);
            fprintf(fid, 'const int32_t %s;\n', bd);

       end
       if(j >= 5)  % a coefficients
            if(j == 5) 
                colstr = '1'; 
            end
            if(j == 6) 
                colstr = '2'; 
            end
            an = ['a' colstr num2str(i) 'n' ' = ' num2str(n(i,j))];
            ad = ['a' colstr num2str(i) 'd' ' = ' num2str(d(i,j))];
            fprintf(fid, 'const int32_t %s;\n', an);
            fprintf(fid, 'const int32_t %s;\n', ad);
       end
   end
end

% write the end of the file
fprintf(fid, '\n\n\n#endif');
fclose(fid);
$\endgroup$

migrated from stackoverflow.com Mar 2 '12 at 0:24

This question came from our site for professional and enthusiast programmers.

  • 4
    $\begingroup$ Higher-order IIR filters like this one are usually implemented using second-order sections; you get the filter you want by cascading multiple second-order stages (with a single first-order stage if the desired order is odd). It is typically a more robust implementation than directly implementing the higher-order filter. $\endgroup$ – Jason R Mar 2 '12 at 1:42
  • 3
    $\begingroup$ If you don't do what @JasonR suggests you will have very large word sizes. Filters like this can fail single precision floating point when implemented with a basic IIR structure, you need the SOS. $\endgroup$ – Christopher Felton Mar 2 '12 at 2:17
  • $\begingroup$ @JasonR: Thank you for suggesting this. I have updated by answer above. $\endgroup$ – Nicholas Kinar Mar 2 '12 at 2:36
  • $\begingroup$ @ChristopherFelton: Thank you for helping to reinforce this. $\endgroup$ – Nicholas Kinar Mar 2 '12 at 2:37
  • $\begingroup$ Yes with your new SOS matrix you can create 3 filters from the site. Or you can use my code here. It will work the same as the website. I will gladly update the script to except SOS matrix. $\endgroup$ – Christopher Felton Mar 2 '12 at 3:19
5
$\begingroup$

As discussed it is best to use the sum of sections, that is break the higher-order filter into a cascaded 2nd order filters. The updated question has the SOS matrix. Using this code and an example here the Python object can be used to generate the individual sections.

In matlab

save SOS

In Python

import shutil
import numpy
from scipy.io import loadmat
from siir import SIIR

matfile = loadmat('SOS.mat')  
SOS = matfile['SOS']
b = numpy.zeros((3,3))
a = numpy.zeros((3,3))
section = [None for ii in range(3)]
for ii in xrange(3):
    b[ii] = SOS[ii,0:3]
    a[ii] = SOS[ii,3:6]

    section[ii] = SIIR(b=b[ii], a=a[ii], W=(24,0))
    section[ii].Convert()  # Create the Verilog for the section
    shutil.copyfile('siir_hdl.v', 'iir_sos_section%d.v'%(ii))

Additional information on fixed-point can be found here

$\endgroup$
  • $\begingroup$ Thank you very much for all of the insightful links, and for the Python code; I hope that your answer (and the other answers posted here) serve as good references for many others. I wish that I could mark all of the responses here as the accepted one. $\endgroup$ – Nicholas Kinar Mar 2 '12 at 14:38
  • 1
    $\begingroup$ If you have any issues let me know and I will update/fix the code if it doesn't work for you. I will modify it (relatively soon, doh) to directly accept a SOS matrix. $\endgroup$ – Christopher Felton Mar 2 '12 at 15:52
  • 1
    $\begingroup$ I tried to implement my own version from your example. On my system, I had to use "from numpy import zeros" and change loatmat to loadmat(). Is the SOS matrix given by Matlab (mathworks.com/help/toolbox/signal/ref/ss2sos.html) in the same format as expected? I receive the following error when trying to access the SOS matrix: "TypeError: unhashable type" when the interpreter reaches the line "b[ii] = SOS[0:3, ii]" $\endgroup$ – Nicholas Kinar Mar 2 '12 at 16:55
  • 1
    $\begingroup$ That would depend on the format of the SOS.mat file. If you simply >>> print(matfile) it will show you the keys in the loaded .mat file. The scipy.io.loadmat always loads as a dictionary (BOMK). $\endgroup$ – Christopher Felton Mar 2 '12 at 18:05
  • 1
    $\begingroup$ Yes, that is correct, the output of 0 is the input to 1 and so forth. A little thought needs to be put in to the word width. The default is two use 24 bit fraction (0 integer, 23 fraction, 1 sign). I believe you originally wanted to use a smaller word width. $\endgroup$ – Christopher Felton Mar 6 '12 at 2:33
10
$\begingroup$

The 'fractional bits' are the number of bits in a bus that you've dedicated to represent the fractional portion of a number (eg, the .75 in 3.75).

Say you've a digital bus 4 bits wide, what number does 1001 represent? It could mean '9' if you treat it as a positive integer (2^3 + 2^0 = 8 + 1 = 9). Or it could mean -7 in two's complement notation: (-2^3 + 2^0 = -8 + 1 = -7).

What about numbers with some fractions in them, ie 'real' numbers? Real numbers can be represented in hardware as "fixed-point" or "floating point". It looks like those filter generators use fixed point.

Back to our 4 bit bus (1001). Let's introduce a binary point so we get 1.001. What this means is that were now using the bits on the RHS of the point to build integers, and the bits on the LHS to build a fraction. The number represented by a digital bus set to 1.001 is 1.125 (1*2^0 + 0*2^-1 + 0*2^-2 + 1*2^-3 = 1 + 0.125 = 1.125). In this case, out of the 4 bits in the bus, we are using 3 of them to represent the fractional portion of a number. Or, we have 3 fractional bits.

So if you've a list of real numbers like you have above, you now have to decide how many fractional bits you want to represent them. And here's the trade-off: the more fractional bits you use, the closer you can represent the number you want, but the bigger your circuit will need to be. And what's more, the fewer fractional bits you use, the further the actual frequency response of the filter will deviate from the one you designed at the start!

And to make matters worse, you're looking to build an Infinite Impulse Response (IIR) filter. These can actually go unstable if you don't have enough fractional and integer bits!

$\endgroup$
  • $\begingroup$ Thanks for providing this insightful answer. I've been trying to run the code generator using the small b coefficients above, and I still receive some errors. Could you suggest something that I could do to properly run the generator? I will update the answer above to show what I've done. $\endgroup$ – Nicholas Kinar Mar 1 '12 at 14:45
10
$\begingroup$

So Marty has well taken care of the bits question. Onto the filter itself, I think you are likely getting a warning or complaint from matlab about poorly scaled coefficients? When I plot the filter, from scipy not matlab but it's likely very similar.

Response

Which is 100 dB down at the passband! So, you may want to make sure you want a smaller order filter, which will help with you implementation anyway. When I get to a 6th order filter I stop getting complaints about bad coefficients. Maybe try reducing the order and see if it still meets your requirements.

$\endgroup$
  • $\begingroup$ Thanks for suggesting this! I think that a 6th order filter would work just as well. Using matlab's fvtool, I think that the response is good for my application. I've now updated my response above. However, something is still going wrong with the Verilog HDL code generator (spiral.net/hardware/filter.html). Perhaps it wants the [b,a] in another format. In addition, +1 for the use of SciPy. $\endgroup$ – Nicholas Kinar Mar 1 '12 at 16:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.