¿Cómo reconstruir datos de fase a partir de datos solo de magnitud en una función de dominio de frecuencia?

Question

Estoy intentando hacer un reverse-engineering de un software llamado LEAP de la ahora desaparecida compañía LinearX. Este software puede realizar una operación llamada "Transformada de Fase Mínima". Esto toma datos reales/magnitud solamente en una curva de voltaje vs frecuencia que se supone es la respuesta de frecuencia de una función de transferencia, y produce la curva de datos complejos completa, agregando los datos de fase. Suponiendo que el sistema es de fase mínima.

Las curvas están en el dominio de frecuencia, con escala/espaciado logarítmico para el eje de frecuencia, el dominio de frecuencia es el rango de frecuencia de audio, de 10Hz a 40000 Hz. Típicamente, las curvas tienen alrededor de 400 puntos de datos.

Las curvas se almacenan y procesan como tablas de valores. Por ejemplo, un filtro de paso de banda ajustado a 1KHz:

Frecuencia   Voltaje 

10 Hz     -40 dBV
10.21 Hz  -39.82 dBV
...
1 KHz      0 dBV
...
40 KHz      -32.04 dBV

El objetivo es tomar una curva como esta:

Curva de voltaje sin datos de fase.

y producir esto:

Curva de voltaje con datos de fase.

Parece que este software utiliza la transformada de Hilbert. Como se indica en Wikipedia:

$$\arg\left[H(j\omega)\right] = -\mathcal{H}\{\log|H(j\omega)|\}$$

Necesita conocer las pendientes/derivadas en ambos extremos de la curva, el extremo de baja frecuencia y el extremo de alta frecuencia en dB/oct. El programa le pide al usuario los dos números de pendiente o intenta estimarlos automáticamente.

Estoy intentando reproducir esta operación con el lenguaje C# en un simple programa de escritorio, o con una simple hoja de cálculo de Excel, pero no he tenido éxito en absoluto.

Mi pregunta es, ¿cuál podría ser el algoritmo para realizar esta operación? Estoy teniendo dificultades para intentar aplicar la transformada de Hilbert. No sé cómo realizar una integral definida impropia desde 0Hz hasta Hz infinito con datos discretos de rango finito.

Parece que la transformada tiene una notación alternativa, que es la convolución de \$f(\omega)\$ y \$g(\omega) = 1/(\pi \omega)\$. Nuevamente, no tengo idea de cómo hacerlo programáticamente. Tal vez podría hacer una convolución discreta, pero ¿qué rango o ventana debo darle a \$g(\omega)\$?

¿Cómo puedo hacer esto? ¿Cómo puedo encontrar la curva de fase a partir de la curva de magnitud, con puntos de datos discretos en el dominio de frecuencia?

Answer 1

Puedes hacerlo con un FFT haciendo que el Cepstrum sea causal (que, hay que admitirlo, es un buen bocado).

1. Tomar el logaritmo natural de la magnitud
2. Hacer una FFT inversa (que es el Cepstrum)
3. Hacerlo causal: duplicar los valores t>0, establecer los valores en t<0 a cero y dejar t=0 y t=N/2 tal como están
4. Hacer una FFT directa. La parte real es el logaritmo de la magnitud, la parte imaginaria es la fase mínima (en radianes).
5. Tomar la exponencial para obtener la función de transferencia en unidades lineales

Hay algunas alternativas a esto:

1. Filtrado en el dominio de frecuencia: la transformada de Hilbert es simplemente un desplazamiento de fase de 90 grados, es decir, la función de transferencia es \$H(\omega) = j*sign(\omega)\$ , donde $sign$ es la función signum
2. Filtrado en el dominio del tiempo. Crear la función de transferencia del paso 1 en una cuadrícula FFT adecuada y hacer una FFT inversa para obtener una respuesta al impulso. Rotar y aplicar una ventana y luego aplicar como filtro FIR.

Hacer esto en una computadora añade la complejidad de que tienes señales discretas en ambos dominios que son periódicas con la frecuencia de muestreo y la longitud del FFT respectivamente y que la multiplicación en el dominio de frecuencia implementa la convolución circular, no la lineal.

Recomendación de lectura sobre la Transformada de Hilbert discreta: http://andrewduncan.net/air/