Eu quero fazer a modelagem de ruído em um aplicativo de 100kHz e 16 bits, para mudar todo o ruído de quantização para a faixa de 25khz-50kHz, com um mínimo de ruído na banda de DC-25kHz.

Configurei o mathematica para criar um kernel de filtro de erro de 31 amostras por meio do aprendizado por reforço que funciona bem: Após um pouco de aprendizado, posso obter ~ 16dB de aprimoramento de ruído de alta frequência para aproximadamente a mesma quantidade de redução na banda de baixa frequência (o linha central é o nível de ruído pontilhado sem forma). Isso está de acordo com o teorema de modelagem de ruído "Gerzon-Craven".

Agora, para o meu problema:

Não consigo modelar o ruído ainda mais depois de um aprendizado extenso, embora o teorema de Gerzon-Craven não o proíba. Por exemplo, deve ser possível obter redução de 40 dB na banda baixa e aprimoramento de 40 dB na banda alta.

Então, existe outro limite fundamental em que estou correndo?

Tentei investigar os teoremas de ruído / amostragem / informação de Shannon, mas, depois de algum tempo, só consegui derivar um único limite: o teorema de Gerzon-Craven, que parece ser um resultado direto do teorema de Shannon.

Qualquer ajuda é apreciada.

EDIT: mais informações

Primeiro, o kernel do filtro que produz a modelagem de ruído acima. Observe que a amostra mais recente está no lado direito . Os valores numéricos do BarChart arredondados para 0,01: {-0,16, 0,51, -0,74, 0,52, -0,04, -0,25, 0,22, -0,11, -0,02, 0,31, -0,56, 0,45, -0,13, 0,04, -0,14, 0,12, -0,06, 0,19, -0,22, -0,15, 0,4, 0,01, -0,41, -0,1, 0,84, -0,42, -0,81, 0,91, 0,75, -2,37, 2,29} (Não é exatamente o caractere de barra, mas produz uma curva semelhante )

Outra observação sobre a implementação do feedback de erro:

Eu tentei duas implementações diferentes do feedback de erro. Primeiro, comparei a amostra de saída arredondada com o valor desejado e usei esse desvio como erro. Segundo, comparei a amostra de saída arredondada com a (entrada + feedback de erro). Embora ambos os métodos produzam kernels bastante diferentes, ambos parecem se nivelar com a mesma intensidade de modelagem de ruído. Os dados publicados aqui usam a segunda implementação.

Aqui está o código usado para calcular as amostras de ondas digitalizadas. step é o tamanho da etapa para o arredondamento. wave é a forma de onda não digitalizada (normalmente apenas zeros quando nenhum sinal é aplicado).

TestWave[kernel_?VectorQ] := 
 Module[{k = kernel, nf, dith, signals, twave, deltas},
  nf = Length@k;
  dith = RandomVariate[TriangularDistribution[{-1, 1}*step], l];
  signals = deltas = Table[0, {l}];
  twave = wave;
  Do[
   twave[[i]] -= k.PadLeft[deltas[[;; i - 1]], nf];
   signals[[i]] = Round[twave[[i]] + dith[[i]], step];
   deltas[[i]] = signals[[i]] - twave[[i]];
   , {i, l}];
  signals
  ]

O método de reforço:

A "pontuação" é calculada observando o espectro de potência do ruído. O objetivo é minimizar a potência do ruído na banda DC-25kHz. Estou não penalizar ruído na banda de alta frequência, de modo que arbitrariamente elevado ruído não haveria pontuação diminuir. Estou introduzindo ruído nos pesos do kernel para aprender. Talvez, portanto, eu esteja em um mínimo local (muito amplo e profundo), mas considero isso extremamente improvável.

Comparação com o design de filtro padrão:

O Mathematica permite gerar filtros iterativamente. Eles podem ter um contraste muito melhor que 36 dB quando a resposta de frequência é plotada; até 80-100 dB. Valores numéricos: {0,024, -0,061, -0,048, 0,38, -0,36, -0,808, 2,09, -0,331, -4,796, 6,142, 3,918, -17,773, 11,245, 30,613, -87,072, 113,676, -87,072, 30,613, 11,245 , -17.773, 3.918, 6.142, -4.796, -0.331, 2.09, -0.808, -0.36, 0.38, -0.048, -0.061, 0.024}

No entanto, ao aplicar aqueles na modelagem de ruído real, eles são (a) fixados no mesmo contraste de ~ 40dB, (b) apresentam desempenho pior do que o filtro aprendido, sem, na verdade, atenuação de ruído.

Pontilhamento básico sem modelagem de ruído

A quantização pontilhada básica sem modelagem de ruído funciona assim:

Figura 1. Diagrama básico do sistema de quantização pontilhada. O ruído é um pontilhamento triangular com média zero e um valor absoluto máximo de 1. O arredondamento é para o número inteiro mais próximo. Erro residual é a diferença entre saída e entrada e é calculado apenas para análise.

$\frac{1}{12}$$\frac{1}{4}$

Com o erro residual aditivo independente, teríamos um modelo mais simples do sistema:

Figura 2. Aproximação da quantização pontilhada básica. Erro residual é ruído branco.

No modelo aproximado, a saída é simplesmente inserida mais erro residual de ruído branco independente.

Pontilhamento com modelagem de ruído

Não consigo ler muito bem o Mathematica, então, em vez do seu sistema, analisarei o sistema de Lipshitz et al. " Modelagem de ruído minimamente audível " J. Audio Eng. Soc., Vol.39, No.11, novembro de 1991:

Figura 3. Lipshitz et al. Diagrama do sistema de 1991 (adaptado da Fig. 1). O filtro (em itálico no texto) inclui um atraso de uma amostra para que possa ser usado como um filtro de feedback de erro. O ruído é pontilhamento triangular.

Se o erro residual é independente dos valores atuais e passados ​​do sinal A, temos um sistema mais simples:

Figura 4. Um modelo aproximado de Lipshitz et al. 1991 sistema. O filtro é o mesmo da Fig. 3 e inclui um atraso de uma amostra. Não é mais usado como um filtro de feedback. Erro residual é ruído branco.

Nesta resposta, trabalharei com o modelo aproximado mais facilmente analisado (Fig. 4). No original Lipshitz et al. Em 1991, o Filter tem uma forma genérica de filtro de resposta ao impulso infinito (IIR) que abrange os filtros IIR e resposta de impulso finito (FIR). A seguir, assumiremos que o Filtro é um filtro FIR, como acredito, com base em meus experimentos com seus coeficientes, que é isso que você tem em seu sistema. A função de transferência do filtro é:

$z^{-1}$

$\ldots, -b_3, -b_2, -b_1$$1, b_1, b_2, b_3, \ldots$$b_0 = 1$horzcatno script Octave abaixo) e, finalmente, a lista é revertida (por flip):

pkg load signal
b = [-0.16, 0.51, -0.74, 0.52, -0.04, -0.25, 0.22, -0.11, -0.02, 0.31, -0.56, 0.45, -0.13, 0.04, -0.14, 0.12, -0.06, 0.19, -0.22, -0.15, 0.4, 0.01, -0.41, -0.1, 0.84, -0.42, -0.81, 0.91, 0.75, -2.37, 2.29];
c = flip(horzcat(-b, 1));
freqz(c)
zplane(c)

O script plota a resposta de frequência de magnitude e os locais zero do filtro de modelagem de ruído completo:

Figura 5. Resposta da frequência de magnitude do filtro de modelagem de ruído completo.

$\times$$\circ$

Penso que o problema de encontrar os coeficientes do filtro pode ser reformulado como o problema de projetar um filtro de fase mínima com um coeficiente inicial de 1. Se houver limitações inerentes à resposta de frequência de tais filtros, essas limitações serão transferidas para limitações equivalentes. na modelagem de ruído que usa esses filtros.

Conversão do projeto de todos os polos para a FIR de fase mínima

Um procedimento para o design de filtros diferentes, mas de muitas maneiras equivalentes, é descrito em Stojanović et al. , "Projeto de filtros digitais recursivos de todos os pólos com base em polinômios ultra-esféricos", Radioengineering, vol. 23, n. 3, setembro de 2014. Eles calculam os coeficientes do denominador da função de transferência de um filtro passa-baixo de todos os pólos IIR. Elas sempre têm um coeficiente denominador inicial de 1 e todos os pólos dentro do círculo unitário, uma exigência de filtros IIR estáveis. Se esses coeficientes forem usados ​​como coeficientes do filtro de modelagem de ruído FIR de fase mínima, eles fornecerão uma resposta de freqüência passa alta invertida em comparação com o filtro IIR passa baixo (os coeficientes denominadores da função de transferência se tornam coeficientes numeradores). Na sua notação {-0.0076120, 0.0960380, -0.5454670, 1.8298040, -3.9884220, 5.8308660, -5.6495140, 3.3816780}, existe um conjunto de coeficientes desse artigo , que pode ser testado para a aplicação de modelagem de ruído, embora não seja exatamente a especificação:

Figura 7. Resposta em frequência de magnitude do filtro FIR usando coeficientes de Stojanović et al. 2014.

Figura 8. Gráfico de pólo zero do filtro FIR usando coeficientes de Stojanović et al. 2014.

A função de transferência de todos os pólos é:

$a$$b$

Para projetar um filtro de todos os pólos e convertê-lo em um filtro FIR de fase mínima, você não poderá usar métodos de design de filtro IIR que iniciem em um filtro de protótipo analógico e mapeie os pólos e zeros para o domínio digital usando a transformação bilinear . Isso inclui cheby1, cheby2e ellipem Octave e SciPy do Python. Esses métodos afastarão os zeros da origem do plano z, para que o filtro não seja do tipo de polo exigido.

Resposta à questão teórica

Se você não se importa com quanto ruído haverá em frequências acima de um quarto da frequência de amostragem, Lipshitz et al. 1991 aborda sua pergunta diretamente:

Para tais funções de ponderação, que chegam a zero sobre parte da banda, não há limite teórico para a redução da potência sonora ponderada obtida no circuito da Fig. 1. Esse seria o caso se, por exemplo, se assumir que o o ouvido tem sensibilidade zero entre, digamos, 20 kHz e a frequência de Nyquist e escolhe a função de ponderação para refletir esse fato.

A Figura 1. mostra um modelador de ruído com uma estrutura de filtro IIR genérica com pólos e zeros, tão diferente da estrutura FIR que você possui no momento, mas o que eles dizem se aplica também a isso, porque uma resposta de impulso do filtro FIR pode ser arbitrariamente próximo da resposta ao impulso de qualquer filtro IIR estável.

Script de oitava para design de filtro

$\nu=0$dip

pkg load signal
N = 14; #number of taps including leading tap with coefficient 1
att = 97.5; #dB attenuation of Dolph-Chebyshev window, must be positive
dip = 2; #spectrum lift-up multiplier, must be above 1
c = chebwin(N, att);
c = conv(c, c);
c /= sum(c);
c(N) += dip*10^(-att/10);
r = roots(c);
j = (abs(r(:)) <= 1);
r = r(j);
c = real(poly(r));
c .*= (-1).^(0:(N-1)); #if this complains, then root finding has probably failed
freqz(c)
zplane(c)
printf('%f, ', flip(-c(2:end))), printf('\n'); #tobalt's format

Começa com uma janela Dolph-Chebyshev como coeficientes, convence-a a dobrar a função de transferência de zeros, adiciona à torneira do meio um número que "eleva" a resposta de frequência (considerando a torneira do meio no tempo zero), de modo que que ele é positivo em todos os lugares, encontra os zeros, remove os zeros que estão fora do círculo unitário, converte os zeros de volta em coeficientes (o coeficiente inicial de polyé sempre 1) e inverte o sinal de cada segundo coeficiente para fazer o filtro passar alto . Os resultados do (uma versão mais antiga, mas quase equivalente) do script parecem promissores:

Figura 9. Resposta da frequência de magnitude do filtro (uma versão mais antiga, mas quase equivalente) do script acima.

Figura 10. Gráfico de pólo zero do filtro (uma versão mais antiga, mas quase equivalente) do script acima.

Os coeficientes de (uma versão mais antiga, mas quase equivalente a) O script acima em sua notação: {0.357662, -2.588396, 9.931419, -26.205448, 52.450624, -83.531276, 108.508775, -116.272581, 102.875781, -74.473956, 43.140431, -19.131434, 5.923468}. Os números são grandes, o que pode levar a problemas numéricos.

Implementação de oitava da modelagem de ruído

Finalmente, eu fiz minha própria implementação de modelagem de ruído no Octave e não tenho problemas como você. Com base em nossa discussão nos comentários, acho que a limitação em sua implementação foi que o espectro de ruído foi avaliado usando uma janela retangular, chamada "sem janelas", que espalhou o espectro de alta frequência nas baixas frequências.

pkg load signal
N = length(c);
M = 16384; #signal length
input = zeros(M, 1);#sin(0.01*(1:M))*127;
er = zeros(M, 1);
output = zeros(M, 1);
for i = 1:M
  A = input(i) + er(i);
  output(i) = round(A + rand() - rand());
  for j = 2:N
    if (i + j - 1 <= M)
      er(i + j - 1) += (output(i) - A)*c(j);
    endif
  endfor
endfor
pwelch(output, max(nuttallwin(1024), 0), 'semilogy');

Figura 11. Análise espectral de ruído de quantização da implementação de oitava acima da Oitava, para modelagem de ruído para sinal de entrada zero constante. Eixo horizontal: frequência normalizada. Preto: sem ruído ( c = [1];), vermelho: seu filtro original, azul: o filtro da seção "Script de oitava para o design do filtro".

Figura 12. Saída no domínio do tempo da implementação Octave acima da modelagem de ruído para o sinal de entrada zero constante. Eixo horizontal: número da amostra, eixo vertical: valor da amostra. Vermelho: seu filtro original, azul: o filtro da seção "Script de oitava para o design do filtro".

O filtro de modelagem de ruído mais extremo (azul) resulta em valores de amostra quantificados de saída muito grandes, mesmo para entrada zero.