Diseñando el filtro Butterworth en Matlab y obteniendo los coeficientes del filtro [ab] como enteros para el generador de código Verilog HDL en línea

15

Diseñé un filtro Butterworth de paso bajo muy simple usando Matlab. El siguiente fragmento de código demuestra lo que he hecho.

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

En [b, a] se almacenan los coeficientes del filtro. Me gustaría obtener [b, a] como enteros para poder usar un generador de código HDL en línea para generar código en Verilog.

Los valores de Matlab [b, a] parecen ser demasiado pequeños para usar con el generador de código en línea (el script Perl del lado del servidor se niega a generar código con los coeficientes), y me pregunto si sería posible obtener [b, a] en una forma que se puede utilizar como una entrada adecuada.

Los coeficientes a que obtengo en Matlab son:

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

Los coeficientes b que obtengo en Matlab son:

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

Usando el generador en línea, me gustaría diseñar un filtro con un ancho de bit de 12 bits y un formulario de filtro I o II. No sé qué se entiende por "bits fraccionales" en el enlace anterior.

Al ejecutar el generador de código (http://www.spiral.net/hardware/filter.html) con los coeficientes [b, a] enumerados anteriormente, con bits fraccionarios establecidos en 20 y un ancho de bits de 12, recibo el siguiente error de ejecución :

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. 

¿Cómo puedo cambiar mi diseño para que este error no ocurra?

ACTUALIZACIÓN: Al usar Matlab para generar un filtro Butterworth de sexto orden, obtengo los siguientes coeficientes:

Para:

1.0000
-5.4914
12.5848
-15.4051
10.6225
-3.9118
0.6010 

para b:

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

Al ejecutar el generador de código en línea (http://www.spiral.net/hardware/filter.html), ahora recibo el siguiente error (con bits fraccionales como 8 y ancho de bits de 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.

¿Quizás los coeficientes b son demasiado pequeños, o quizás el generador de código (http://www.spiral.net/hardware/filter.html) quiere el [b, a] en otro formato?

ACTUALIZAR:

Quizás lo que necesito hacer es escalar los coeficientes [b, a] por el número de bits fraccionarios para obtener los coeficientes como enteros.

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

Sin embargo, sigo pensando que los coeficientes b son extremadamente pequeños. ¿Qué estoy haciendo mal aquí?

¿Quizás otro tipo de filtro (o método de diseño de filtro) sería más adecuado? ¿Alguien podría hacer una sugerencia?

ACTUALIZACIÓN: Según lo sugerido por Jason R y Christopher Felton en los comentarios a continuación, un filtro SOS sería más adecuado. Ahora he escrito un código de Matlab para obtener un filtro SOS.

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);

La matriz SOS que obtengo es:

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

¿Es posible seguir utilizando la herramienta de generación de código Verilog (http://www.spiral.net/hardware/filter.html) para implementar este filtro SOS, o debería simplemente escribir el Verilog a mano? ¿Hay una buena referencia disponible?

Me pregunto si un filtro FIR sería mejor para usar en esta situación.

MÁS INFORMACIÓN: Los filtros IIR recursivos se pueden implementar usando matemáticas enteras expresando coeficientes como fracciones. (Consulte el excelente libro de procesamiento de señales DSP de Smith para obtener más detalles: http://www.dspguide.com/ch19/5.htm )

El siguiente programa de Matlab convierte los coeficientes de filtro de Butterworth en partes fraccionarias utilizando la función Matlab rat (). Luego, como se menciona en los comentarios, las secciones de segundo orden se pueden usar para implementar numéricamente el filtro (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);
Nicholas Kinar
fuente
44
Los filtros IIR de orden superior como este generalmente se implementan mediante secciones de segundo orden ; obtienes el filtro que deseas conectando en cascada múltiples etapas de segundo orden (con una sola etapa de primer orden si el orden deseado es impar). Por lo general, es una implementación más sólida que la implementación directa del filtro de orden superior.
Jason R
3
Si no hace lo que @JasonR sugiere, tendrá un tamaño de palabra muy grande. Los filtros como este pueden fallar con coma flotante de precisión simple cuando se implementan con una estructura básica IIR, necesita el SOS.
Christopher Felton
@JasonR: Gracias por sugerir esto. He actualizado por respuesta anterior.
Nicholas Kinar
@ChristopherFelton: Gracias por ayudar a reforzar esto.
Nicholas Kinar
Sí, con su nueva matriz SOS puede crear 3 filtros desde el sitio. O puedes usar mi código aquí . Funcionará igual que el sitio web. Con mucho gusto actualizaré el script a excepción de la matriz SOS.
Christopher Felton

Respuestas:

5

Como se discutió, es mejor usar la suma de las secciones, es decir, dividir el filtro de orden superior en filtros de segundo orden en cascada. La pregunta actualizada tiene la matriz SOS. Usando este código y un ejemplo aquí, el objeto Python se puede usar para generar las secciones individuales.

En matlab

save SOS

En 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))

Puede encontrar información adicional sobre punto fijo aquí

Christopher Felton
fuente
Muchas gracias por todos los enlaces reveladores y por el código Python; Espero que su respuesta (y las otras respuestas publicadas aquí) sirvan como buenas referencias para muchos otros. Desearía poder marcar todas las respuestas aquí como aceptadas.
Nicholas Kinar
1
Si tiene algún problema, hágamelo saber y actualizaré / arreglaré el código si no funciona para usted. Lo modificaré (relativamente pronto, doh) para aceptar directamente una matriz SOS.
Christopher Felton
1
Traté de implementar mi propia versión a partir de su ejemplo. En mi sistema, tuve que usar "from numpy import ceros" y cambiar loatmat a loadmat (). ¿Está la matriz SOS dada por Matlab ( mathworks.com/help/toolbox/signal/ref/ss2sos.html ) en el mismo formato que se esperaba? Recibo el siguiente error al intentar acceder a la matriz SOS: "TypeError: tipo no compartible" cuando el intérprete llega a la línea "b [ii] = SOS [0: 3, ii]"
Nicholas Kinar
1
Eso dependería del formato del archivo SOS.mat. Si simplemente >>> imprime (matfile) le mostrará las claves en el archivo .mat cargado. Scipy.io.loadmat siempre se carga como un diccionario (BOMK).
Christopher Felton
1
Sí, eso es correcto, la salida de 0 es la entrada a 1 y así sucesivamente. Hay que pensar un poco en el ancho de la palabra. El valor predeterminado es dos fracciones de 24 bits (0 enteros, 23 fracciones, 1 signo). Creo que originalmente querías usar un ancho de palabra más pequeño.
Christopher Felton
10

Los 'bits fraccionarios' son el número de bits en un bus que ha dedicado a representar la parte fraccionaria de un número (por ejemplo, el .75 en 3.75).

Digamos que tiene un bus digital de 4 bits de ancho, ¿qué número 1001representa? Podría significar '9' si lo trata como un entero positivo (2 ^ 3 + 2 ^ 0 = 8 + 1 = 9). O podría significar -7 en notación de complemento a dos: (-2 ^ 3 + 2 ^ 0 = -8 + 1 = -7).

¿Qué pasa con los números con algunas fracciones en ellos, es decir, números 'reales'? Los números reales se pueden representar en hardware como "punto fijo" o "punto flotante". Parece que esos generadores de filtro usan punto fijo.

De vuelta a nuestro bus de 4 bits ( 1001). Vamos a introducir un punto binario para obtener 1.001. Lo que esto significa es que ahora estábamos usando los bits en el RHS del punto para construir enteros, y los bits en el LHS para construir una fracción. El número representado por un bus digital configurado 1.001es 1.125 ( 1* 2 ^ 0 + 0* 2 ^ -1 + 0* 2 ^ -2 + 1* 2 ^ -3 = 1 + 0.125 = 1.125). En este caso, de los 4 bits en el bus, estamos usando 3 de ellos para representar la porción fraccional de un número. O tenemos 3 bits fraccionarios.

Entonces, si tienes una lista de números reales como la que tienes arriba, ahora tienes que decidir cuántos bits fraccionarios quieres representar. Y aquí está la compensación: cuantos más bits fraccionales use, más cerca podrá representar el número que desea, pero más grande deberá ser su circuito. ¡Y lo que es más, cuantos menos bits fraccionarios use, más se desviará la respuesta de frecuencia real del filtro del que diseñó al principio!

Y para empeorar las cosas, está buscando construir un filtro de respuesta de impulso infinito (IIR). En realidad, estos pueden volverse inestables si no tienes suficientes bits fraccionales y enteros.

Marty
fuente
Gracias por proporcionar esta respuesta perspicaz. He estado tratando de ejecutar el generador de código usando los coeficientes b pequeños anteriores, y todavía recibo algunos errores. ¿Podría sugerirme algo que podría hacer para ejecutar correctamente el generador? Actualizaré la respuesta anterior para mostrar lo que he hecho.
10

Así que Marty se ha ocupado bien de la cuestión de los bits. En el filtro en sí, ¿creo que es probable que reciba una advertencia o una queja de matlab sobre coeficientes escalados? Cuando trazo el filtro, desde scipy no matlab, pero probablemente sea muy similar.

Respuesta

¡Que está 100 dB abajo en la banda de paso! Por lo tanto, es posible que desee asegurarse de que desea un filtro de orden más pequeño, que de todos modos lo ayudará con su implementación. Cuando llego a un filtro de sexto orden, dejo de recibir quejas sobre malos coeficientes. Tal vez intente reducir el pedido y vea si aún cumple con sus requisitos.

macduff
fuente
¡Gracias por sugerir esto! Creo que un filtro de sexto orden funcionaría igual de bien. Usando fvtool de matlab, creo que la respuesta es buena para mi aplicación. Ahora he actualizado mi respuesta anterior. Sin embargo, algo sigue yendo mal con el generador de código Verilog HDL ( spiral.net/hardware/filter.html ). Quizás quiera la [b, a] en otro formato. Además, +1 por el uso de SciPy.