Projetando o filtro Butterworth no Matlab e obtendo coeficientes do filtro [ab] como números inteiros para o gerador de código Verilog HDL online

15

Projetei um filtro Butterworth low-pass muito simples usando o Matlab. O seguinte trecho de código demonstra o que eu fiz.

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

Em [b, a] são armazenados os coeficientes do filtro. Eu gostaria de obter [b, a] como números inteiros para que eu possa usar um gerador de código HDL online para gerar código no Verilog.

Os valores do Matlab [b, a] parecem muito pequenos para serem usados ​​com o gerador de código on-line (o script Perl do lado do servidor se recusa a gerar código com os coeficientes), e estou pensando se seria possível obter [b, a] de uma forma que possa ser usada como uma entrada adequada.

Os coeficientes a que recebo no Matlab são:

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

Os coeficientes b que recebo no Matlab são:

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 o gerador on-line, eu gostaria de criar um filtro com uma largura de bits de 12 bits e um formulário de filtro I ou II. Não sei o que se entende por "bits fracionários" no link acima.

Executando o gerador de código (http://www.spiral.net/hardware/filter.html) com os coeficientes [b, a] listados acima, com bits fracionários definidos em 20 e uma largura de bit de 12, recebo o seguinte erro de execução :

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. 

Como posso alterar meu design para que esse erro não ocorra?

ATUALIZAR: Usando o Matlab para gerar um filtro Butterworth de 6ª ordem, obtenho os seguintes 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

Executando o gerador de código online (http://www.spiral.net/hardware/filter.html), agora recebo o seguinte erro (com bits fracionários como 8 e largura de bit 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.

Talvez os coeficientes b sejam muito pequenos, ou talvez o gerador de código (http://www.spiral.net/hardware/filter.html) queira o [b, a] em outro formato?

ATUALIZAR:

Talvez o que eu precise fazer seja dimensionar os coeficientes [b, a] pelo número de bits fracionários para obter os coeficientes como números inteiros.

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

No entanto, ainda acho que os coeficientes b são extremamente pequenos. O que eu estou fazendo errado aqui?

Talvez outro tipo de filtro (ou método de design de filtro) seja mais adequado? Alguém poderia fazer uma sugestão?

ATUALIZAÇÃO: Conforme sugerido por Jason R e Christopher Felton nos comentários abaixo, um filtro SOS seria mais adequado. Agora eu escrevi um código Matlab para obter um 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);

A matriz SOS que recebo é:

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

Ainda é possível usar a ferramenta de geração de código Verilog (http://www.spiral.net/hardware/filter.html) para implementar esse filtro SOS, ou devo simplesmente escrever o Verilog à mão? Existe uma boa referência disponível?

Gostaria de saber se um filtro FIR seria melhor usar nesta situação.

Além disso, os filtros IIR recursivos podem ser implementados usando matemática inteira, expressando coeficientes como frações. (Consulte o excelente livro de processamento de sinal DSP da Smith para obter mais detalhes: http://www.dspguide.com/ch19/5.htm )

O programa Matlab a seguir converte os coeficientes do filtro Butterworth em partes fracionárias usando a função rat () do Matlab. Como mencionado nos comentários, as seções de segunda ordem podem ser usadas para implementar numericamente o 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
fonte
4
Filtros IIR de ordem superior como essa geralmente são implementados usando seções de segunda ordem ; você obtém o filtro desejado em cascata de vários estágios de segunda ordem (com um único estágio de primeira ordem, se a ordem desejada for ímpar). Geralmente, é uma implementação mais robusta do que a implementação direta do filtro de ordem superior.
Jason R
3
Se você não fizer o que @JasonR sugere, terá tamanhos de palavras muito grandes. Filtros como esse podem falhar em um ponto flutuante de precisão única quando implementados com uma estrutura IIR básica, você precisa do SOS.
9788 Christopher
@JasonR: Obrigado por sugerir isso. Eu atualizei pela resposta acima.
Nicholas Kinar
@ChristopherFelton: Obrigado por ajudar a reforçar isso.
Nicholas Kinar
Sim, com sua nova matriz SOS, você pode criar 3 filtros a partir do site. Ou você pode usar meu código aqui . Funcionará da mesma forma que o site. Terei prazer em atualizar o script, exceto a matriz SOS.
9309 Christopher

Respostas:

5

Conforme discutido, é melhor usar a soma das seções, que é dividir o filtro de ordem superior em filtros de 2ª ordem em cascata. A pergunta atualizada tem a matriz SOS. Usando esse código e um exemplo aqui, o objeto Python pode ser usado para gerar as seções individuais.

No matlab

save SOS

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

Informações adicionais sobre ponto fixo podem ser encontradas aqui

Christopher Felton
fonte
Muito obrigado por todos os links perspicazes e pelo código Python; Espero que sua resposta (e as outras respostas postadas aqui) sirva como boas referências para muitas outras. Eu gostaria de poder marcar todas as respostas aqui como as aceitas.
Nicholas Kinar
11
Se você tiver algum problema, avise-me e atualizarei / corrigirei o código se ele não funcionar para você. Vou modificá-lo (relativamente cedo, doh) para aceitar diretamente uma matriz SOS.
9309 Christopher
11
Eu tentei implementar minha própria versão do seu exemplo. No meu sistema, tive que usar "from numpy import zeros" e alterar loatmat para loadmat (). A matriz SOS é fornecida pelo Matlab ( mathworks.com/help/toolbox/signal/ref/ss2sos.html ) no mesmo formato que o esperado? Eu recebo o seguinte erro ao tentar acessar a matriz SOS: "TypeError: unhashable type" quando o intérprete atinge a linha "b [ii] = SOS [0: 3, ii]"
Nicholas Kinar
11
Isso dependeria do formato do arquivo SOS.mat. Se você simplesmente >>> imprimir (matfile), ele mostrará as chaves no arquivo .mat carregado. O scipy.io.loadmat sempre é carregado como um dicionário (BOMK).
Christopher Felton
11
Sim, isso está correto, a saída de 0 é a entrada para 1 e assim por diante. Um pouco de reflexão precisa ser colocado na palavra largura. O padrão é usar fração de 24 bits (0 inteiro, fração 23, 1 sinal). Eu acredito que você originalmente queria usar uma largura menor de palavra.
9309 Christopher Felton
10

Os 'bits fracionários' são o número de bits em um barramento que você dedicou para representar a parte fracionária de um número (por exemplo, 0,75 em 3,75).

Digamos que você tenha um barramento digital de 4 bits de largura, que número faz 1001 representa? Pode significar '9' se você o tratar como um número inteiro positivo (2 ^ 3 + 2 ^ 0 = 8 + 1 = 9). Ou pode significar -7 na notação do complemento de dois: (-2 ^ 3 + 2 ^ 0 = -8 + 1 = -7).

E os números com algumas frações, ou seja, números 'reais'? Números reais podem ser representados no hardware como "ponto fixo" ou "ponto flutuante". Parece que esses geradores de filtro usam ponto fixo.

Voltar ao nosso barramento de 4 bits ( 1001). Vamos introduzir um ponto binário para que possamos entender 1.001. O que isto significa é que agora estavam usando os bits no RHS do ponto para construir números inteiros e os bits no LHS para construir uma fração. O número representado por um barramento digital definido como 1.0011.125 ( 1* 2 ^ 0 + 0* 2 ^ -1 + 0* 2 ^ -2 + 1* 2 ^ -3 = 1 + 0,125 = 1,125). Nesse caso, dos 4 bits no barramento, estamos usando 3 deles para representar a parte fracionária de um número. Ou, temos 3 bits fracionários.

Portanto, se você tem uma lista de números reais como você tem acima, agora precisa decidir quantos bits fracionários deseja representá-los. E aqui está a troca: quanto mais bits fracionários você usar, mais perto poderá representar o número que deseja, mas maior será o seu circuito. Além disso, quanto menos bits fracionários você usar, mais a resposta de frequência real do filtro se desviará da que você projetou no início!

E para piorar as coisas, você deseja criar um filtro de resposta ao impulso infinito (IIR). Isso pode ficar instável se você não tiver bits fracionários e inteiros suficientes!

Marty
fonte
Obrigado por fornecer esta resposta perspicaz. Eu tenho tentado executar o gerador de código usando os pequenos coeficientes b acima, e ainda recebo alguns erros. Você poderia sugerir algo que eu poderia fazer para executar corretamente o gerador? Vou atualizar a resposta acima para mostrar o que fiz.
10

Então Marty cuidou bem da questão dos bits. No próprio filtro, acho que você provavelmente está recebendo um aviso ou reclamação do matlab sobre coeficientes com escala insuficiente? Quando plotei o filtro, do scipy, não do matlab, mas provavelmente é muito parecido.

Resposta

Que é de 100 dB na faixa de passagem! Portanto, convém ter um filtro de pedidos menor, o que ajudará na sua implementação de qualquer maneira. Quando chego a um filtro de 6ª ordem, paro de receber reclamações sobre coeficientes ruins. Talvez tente reduzir o pedido e ver se ele ainda atende aos seus requisitos.

macduff
fonte
Obrigado por sugerir isso! Eu acho que um filtro de 6ª ordem também funcionaria. Usando o fvtool do matlab, acho que a resposta é boa para minha aplicação. Atualizei minha resposta acima. No entanto, algo ainda está errado com o gerador de código HDL da Verilog ( spiral.net/hardware/filter.html ). Talvez ele queira o [b, a] em outro formato. Além disso, +1 para o uso do SciPy.