Quiero hacer un modelado de ruido en una aplicación de 16 bits a 100 kHz, para cambiar todo el ruido de cuantificación a la banda de 25 kHz a 50 kHz, con un ruido mínimo en la banda de CC a 25 kHz.

Configuré Mathica para crear un núcleo de filtro de error de 31 muestras a través del aprendizaje de refuerzo que funciona bien: después de un poco de aprendizaje, puedo obtener una mejora de aproximadamente 16 dB de ruido de alta frecuencia para aproximadamente la misma cantidad de reducción en la banda de baja frecuencia (el la línea central es el nivel de ruido de oscilación sin forma). Esto está en línea con el teorema de formación de ruido "Gerzon-Craven".

Ahora a mi problema:

No puedo dar forma al ruido aún más, incluso después de un aprendizaje extenso, aunque el teorema de Gerzon-Craven no lo prohíbe. Por ejemplo, debería ser posible lograr una reducción de 40 dB en la banda baja y una mejora de 40 dB en la banda alta.

Entonces, ¿hay otro límite fundamental con el que me encuentro?

Traté de analizar los teoremas de ruido / muestreo / información de Shannon, pero después de jugar un rato con eso, solo logré derivar un límite: el teorema de Gerzon-Craven, que parece ser un resultado directo del teorema de Shannon.

Cualquier ayuda es apreciada.

EDITAR: más información

En primer lugar, el núcleo del filtro que produce la formación de ruido anterior. Tenga en cuenta que la muestra más reciente está en el lado derecho . Valores numéricos del gráfico de barras redondeados a .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} (No es exactamente la barra de caracteres pero produce una curva similar )

Otra nota sobre la implementación de comentarios de error:

Intenté dos implementaciones diferentes de la retroalimentación de error. Primero, comparé la muestra de salida redondeada con el valor deseado y usé esta desviación como error. En segundo lugar, comparé la muestra de salida redondeada con la (entrada + retroalimentación de error). Aunque ambos métodos producen núcleos bastante diferentes, ambos parecen nivelarse con aproximadamente la misma intensidad de conformación de ruido. Los datos publicados aquí utilizan la segunda implementación.

Aquí está el código que se utiliza para calcular las muestras de onda digitalizadas. step es el tamaño de paso para redondear. wave es la forma de onda no digitalizada (generalmente solo ceros cuando no se aplica señal).

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
  ]

El método de refuerzo:

El "puntaje" se calcula observando el espectro de potencia de ruido. El objetivo es minimizar la potencia de ruido en la banda DC-25kHz. Estoy no penalizar el ruido en la banda de alta frecuencia, por lo que arbitrariamente alto ruido no habría puntuación disminuirá. Estoy introduciendo ruido en los pesos del núcleo para aprender. Quizás, por lo tanto, estoy en un mínimo local (muy amplio y profundo), pero considero que esto es extremadamente improbable.

Comparación con el diseño de filtro estándar:

Mathematica permite generar filtros de forma iterativa. Estos pueden tener un contraste mucho mejor que 36 dB cuando se traza su respuesta de frecuencia; hasta 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}

Sin embargo, cuando se aplican aquellos en la configuración de ruido real, (a) se sujetan al mismo contraste de ~ 40dB, (b) funcionan peor que el filtro aprendido sin realmente atenuar el ruido.

Dithering básico sin conformación de ruido

La cuantización básica difuminada sin conformación de ruido funciona así:

Figura 1. Diagrama básico del sistema de cuantificación diferida. El ruido es un tramado triangular de media cero con un valor absoluto máximo de 1. El redondeo es al entero más cercano. El error residual es la diferencia entre salida y entrada, y se calcula solo para análisis.

El tramado triangular aumenta la varianza del error residual resultante en un factor de 3 (de a ) pero desacopla la media y la varianza del error de cuantificación neta del valor de la señal de entrada. Eso significa que la señal de error neta no está correlacionada con la entrada, pero los momentos más altos no están desacoplados, por lo que no es un error aleatorio completamente independiente, pero nadie ha determinado que las personas puedan escuchar cualquier dependencia de los momentos más altos en la señal de error neta en el señal de entrada en una aplicación de audio.$\frac{1}{12}$$\frac{1}{4}$

Con un error residual aditivo independiente tendríamos un modelo más simple del sistema:

Figura 2. Aproximación de la cuantización básica difusa. El error residual es ruido blanco.

En el modelo aproximado, la salida es simplemente entrada más error residual de ruido blanco independiente.

Dithering con forma de ruido

No puedo leer Mathematica muy bien, así que en lugar de su sistema, analizaré el sistema de Lipshitz et al. " Forma de ruido mínimamente audible " J. Audio Eng. Soc., Vol.39, No.11, noviembre de 1991:

Figura 3. Lipshitz et al. Diagrama del sistema de 1991 (adaptado de su Fig. 1). El filtro (en cursiva en el texto) incluye una demora de una muestra para que pueda usarse como un filtro de retroalimentación de error. El ruido es dither triangular.

Si el error residual es independiente de los valores actuales y pasados ​​de la señal A, tenemos un sistema más simple:

Figura 4. Un modelo aproximado de Lipshitz et al. Sistema de 1991. El filtro es el mismo que en la Fig. 3 e incluye en él una demora de una muestra. Ya no se usa como filtro de retroalimentación. El error residual es ruido blanco.

En esta respuesta trabajaré con el modelo aproximado más fácil de analizar (Fig. 4). En el original Lipshitz et al. Sistema 1991, Filter tiene una forma de filtro genérico de respuesta de impulso infinito (IIR) que cubre tanto los filtros de respuesta de impulso finito (FIR) como IIR. A continuación, supondremos que Filter es un filtro FIR, ya que creo, según mis experimentos con sus coeficientes, que eso es lo que tiene en su sistema. La función de transferencia de Filter es:

El factor representa un retraso de una muestra. En el modelo aproximado también hay una ruta de suma directa a la salida del error residual. Esto se resume con la salida negada del filtro , formando la función de transferencia del filtro de conformación de ruido completo:$z^{-1}$

Para ir desde sus coeficientes de filtro , que enumera en orden , a los coeficientes polinomiales la función de transferencia de filtro de conformación de ruido completo , el signo de los coeficientes se cambia para tener en cuenta la negación de la salida del filtro en el diagrama del sistema, y ​​el coeficiente se agrega al final (por en el script de octava a continuación), y finalmente se invierte la lista (por ): 1 , b 1 , b 2 , b 3 , … b 0 = $\ldots, -b_3, -b_2, -b_1$$1, b_1, b_2, b_3, \ldots$$b_0 = 1$horzcatflip

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)

El script traza la respuesta de frecuencia de magnitud y las ubicaciones cero del filtro de conformación de ruido completo:

Figura 5. Respuesta de frecuencia de magnitud del filtro de conformación de ruido completo.

Figura 6. Gráfico del plano Z de polos ( ) y ceros ( ) del filtro. Todos los ceros están dentro del círculo unitario, por lo que el filtro de conformación de ruido completo es de fase mínima.$\times$$\circ$

Creo que el problema de encontrar los coeficientes de filtro puede reformularse como el problema de diseñar un filtro de fase mínima con un coeficiente principal de 1. Si existen limitaciones inherentes a la respuesta de frecuencia de dichos filtros, estas limitaciones se transfieren a limitaciones equivalentes en forma de ruido que usa tales filtros.

Conversión del diseño de todos los polos a FIR de fase mínima

Un procedimiento para el diseño de filtros diferentes pero en muchos sentidos equivalentes se describe en Stojanović et al. , "Diseño de filtros digitales recursivos de todos los polos basado en polinomios ultraesféricos", Radioengineering, vol 23, no 3, septiembre de 2014. Calculan los coeficientes del denominador de la función de transferencia de un filtro de paso bajo de todos los polos IIR. Esos siempre tienen un coeficiente de denominador principal de 1 y tienen todos los polos dentro del círculo unitario, un requisito de filtros IIR estables. Si esos coeficientes se usan como coeficientes del filtro de conformación de ruido FIR de fase mínima, darán una respuesta de frecuencia de paso alto invertida en comparación con el filtro IIR de paso bajo (los coeficientes del denominador de la función de transferencia se convierten en coeficientes numeradores). En su notación, un conjunto de coeficientes de ese artículo es {-0.0076120, 0.0960380, -0.5454670, 1.8298040, -3.9884220, 5.8308660, -5.6495140, 3.3816780}, que podría probarse para la aplicación de modelado de ruido, aunque no es exactamente la especificación:

Figura 7. Respuesta de frecuencia de magnitud del filtro FIR utilizando coeficientes de Stojanović et al. 2014.

Figura 8. Gráfico de polo cero del filtro FIR utilizando coeficientes de Stojanović et al. 2014.

La función de transferencia de todos los polos es:

Por lo tanto, se puede diseñar un establo IIR filtro paso-bajo de todos los polos y utilizar los coeficientes coeficientes para obtener un filtro de paso alto FIR de fase mínima con un coeficiente principal de 1.$a$$b$

Para diseñar un filtro de todos los polos y convertirlo en un filtro FIR de fase mínima, no podrá utilizar métodos de diseño de filtros IIR que comiencen desde un filtro prototipo analógico y mapeen los polos y ceros en el dominio digital usando la transformación bilineal . Eso incluye cheby1, cheby2y ellipen Octave and Python's SciPy. Estos métodos darán ceros lejos del origen del plano z para que el filtro no sea del tipo de todos los polos requerido.

Respuesta a la pregunta teórica.

Si no le importa cuánto ruido habrá en frecuencias superiores a la cuarta parte de la frecuencia de muestreo, entonces Lipshitz et al. 1991 aborda su pregunta directamente:

Para tales funciones de ponderación, que van a cero en una parte de la banda, no existe un límite teórico para la reducción de potencia de ruido ponderada obtenible del circuito de la Fig. 1. Este sería el caso si, por ejemplo, se supone que el ear tiene una sensibilidad cero entre, digamos, 20 kHz y la frecuencia de Nyquist, y elige la función de ponderación para reflejar este hecho.

Su Fig 1. muestra un modelador de ruido con una estructura de filtro IIR genérica con polos y ceros, tan diferente a la estructura FIR que tiene en este momento, pero lo que dicen también se aplica a eso, porque una respuesta de impulso del filtro FIR puede ser hecho arbitrariamente cerca de la respuesta al impulso de cualquier filtro IIR estable dado.

Script de octava para diseño de filtro

Aquí hay una secuencia de comandos Octave para el cálculo del coeficiente por otro método que creo que es equivalente al Stojanovici et al. Método 2014 parametrizado como con la elección correcta de mi parámetro.$\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

Comienza con una ventana Dolph-Chebyshev como los coeficientes, lo convoluciona consigo mismo para duplicar los ceros de la función de transferencia, agrega al toque medio un número que "eleva" la respuesta de frecuencia (considerando que el toque medio está en tiempo cero) que es positivo en todas partes, encuentra los ceros, elimina los ceros que están fuera del círculo unitario, convierte los ceros nuevamente en coeficientes (el coeficiente inicial de polysiempre es 1) y cambia el signo de cada segundo coeficiente para hacer que el filtro pase alto . Los resultados de (una versión anterior pero casi equivalente) del script parecen prometedores:

Figura 9. Respuesta de frecuencia de magnitud del filtro de (una versión anterior pero casi equivalente de) el script anterior.

Figura 10. Gráfico de polo cero del filtro de (una versión anterior pero casi equivalente) del script anterior.

Los coeficientes de (una versión más antigua, pero casi equivalente a) la secuencia de comandos en su notación: {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}. Los números son grandes, lo que podría conducir a problemas numéricos.

Implementación de octava de modelado de ruido

Finalmente, hice mi propia implementación de modelado de ruido en Octave y no tengo problemas como tú. En base a nuestra discusión en los comentarios, creo que la limitación en su implementación fue que el espectro de ruido se evaluó utilizando una ventana rectangular, también conocida como "sin ventanas", que derramó el espectro de alta frecuencia a las bajas frecuencias.

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álisis espectral de ruido de cuantificación de la implementación de octava anterior de conformación de ruido para señal de entrada cero constante. Eje horizontal: frecuencia normalizada. Negro: sin forma de ruido ( c = [1];), rojo: su filtro original, azul: el filtro de la sección "Script de octava para el diseño del filtro".

Figura 12. Salida en el dominio del tiempo de la implementación de octava anterior de conformación de ruido para señal de entrada cero constante. Eje horizontal: número de muestra, eje vertical: valor de muestra. Rojo: su filtro original, azul: el filtro de la sección "Script de octava para diseño de filtro".

El filtro de conformación de ruido más extremo (azul) da como resultado valores de muestra de salida cuantificados muy grandes incluso para entrada cero.