4
$\begingroup$

I have a filter/LTI system frequency response in form of list of values in the frequency domain. I want to get the phase curve/data from magnitude data.

Input data can have either linear spaced points or log spaced points, in the frequency domain. Data is like following table:

Freq (Hz)     Gain (dB)
 512.588      -1.01   
 529.430      -1.07   
 546.826      -1.14  
 ...

Curves obtained from input data can be very irregular and noisy, but is minimum phase for sure. In wikipedia, is stated the relationship between gain and phase for a minimum phase system is:

\begin{align} \arg \left[ H(j \omega) \right] = -\mathcal{H} \lbrace \ln \left( |H(j \omega)| \right) \rbrace \end{align}

The goal is to get the $ \arg \left[ H(j \omega) \right] $ function. Assuming system is minimum phase. I know this is possible because I'm trying to replicate this operation from an existing software called LinearX LEAP. What I'm trying to replicate is the algorithm.

Obtaining phase data

FFT based solutions are NOT allowed because it would alter the log spacing of the data points of frequency axis. I tried to develop the solution as far as I could but got stuck.

Hilbert Transform is defined as:

\begin{align} \mathcal{H} \lbrace H(\omega) \rbrace \ \stackrel{\mathrm{def}}{=}\ \widehat{H}(\omega) = \frac{1}{\pi}\int_{-\infty}^{\infty}\frac{H(\tau)}{\omega -\tau}\, d\tau \end{align}

So I have to transform the natural logarithm of the data. Like this:

\begin{align} \arg \left[ H(j \omega) \right] = - \frac{1}{\pi}\int_{0}^{\infty}\frac{ln(H(\tau ))}{\omega -\tau}\, d\tau \end{align}

I removed the absolute value because data is real values only, not complex.

Frequency response data is contained in between a minimum frequency $ \omega_L $ and a maximum frequency $ \omega_R $

Layout of data

In order to realize the integration in the full domain of 0 Hz to Infinity Hz, I have separated the integral in three sections:

  1. Low frequency tail: $ H(\omega) = a\omega^{L1} $ with $ 0 < \omega < \omega_L $
  2. Actual data with $ \omega_L < \omega < \omega_R $
  3. High frequency tail: $ H(\omega) = b\omega^{L2} $ with $ \omega_R < \omega < \infty $

Where $ L1 $ is the low frequency slope, $ L2 $ is the high frequency slope. $ a $ and $ b $ are coefficients for splicing the tails with the data. The actual domain for $ \omega $ is $ \omega_L < \omega < \omega_R $. Because the data domain is the frequency range of interest.

Solving the integral for the Low Frequency tail is done as:

\begin{align} \int_{0}^{\omega_L}\frac{ln(a \tau ^ {L1} )}{\omega - \tau}\, d\tau = {L1} \int_{0}^{\omega_L}\frac{ln( \tau )}{\omega - \tau}\, d\tau \end{align}

I remove the coefficient $ a $ because the transform of a constant is zero.

Using Cauchy Principal Value:

\begin{align} = \lim_{\epsilon \to 0} {L1} \int_{\epsilon}^{\omega_L}\frac{ln( \tau )}{\omega - \tau}\, d\tau \end{align}

I got the anti-derivative using WolframAlpha Online Integral Calculator

\begin{align} = \lim_{\epsilon \to 0} L1 \bigg| - Li_2(\frac{\tau}{\omega}) - ln(\tau) ln(1 - \frac{\tau}{\omega}) \bigg|_{\epsilon}^{\omega_L} \end{align}

\begin{align} = \lim_{\epsilon \to 0} L1 \left( - Li_2(\frac{\omega_L}{\omega}) - ln(\omega_L) ln(1 - \frac{\omega_L}{\omega}) + Li_2(\frac{\epsilon}{\omega}) + ln(\epsilon) ln(1 - \frac{\epsilon}{\omega}) \right) \end{align}

\begin{align} = L1 \left( - Li_2(\frac{\omega_L}{\omega}) - ln(\omega_L) ln(1 - \frac{\omega_L}{\omega}) \right) \end{align}

So far, everything is fine. The argument of the $ Li_2 $ function is lesser than one, because $ \omega_L < \omega $

The second part of the integral is integrated with numerical integration. We have to use the Cauchy Principal Value where $ \tau = \omega $ . No problem there.

Now in the third part of the integral come the problems. This is integrated between $ \omega_R $ and $ \infty $ . I do:

\begin{align} \int_{\omega_R}^{\infty}\frac{ln(b \tau ^ {L2} )}{\omega - \tau}\, d\tau = {L2} \int_{\omega_R}^{\infty}\frac{ln( \tau )}{\omega - \tau}\, d\tau \end{align}

Using Principal Value

\begin{align} = \lim_{c \to \infty} {L2} \int_{\omega_R}^{c}\frac{ln( \tau )}{\omega - \tau}\, d\tau \end{align}

\begin{align} = \lim_{c \to \infty} L2 \bigg| - Li_2(\frac{\tau}{\omega}) - ln(\tau) ln(1 - \frac{\tau}{\omega}) \bigg|_{\omega_R}^{c} \end{align}

\begin{align} = \lim_{c \to \infty} L2 \left( - Li_2(\frac{c}{\omega}) - ln(c) ln(1 - \frac{c}{\omega}) + Li_2(\frac{\omega_R}{\omega}) + ln(\omega_R) ln(1 - \frac{\omega_R}{\omega}) \right) \end{align}

The terms $ - Li_2(\frac{c}{\omega}) $ and $ - ln(c) ln(1 - \frac{c}{\omega}) $ seems to tend to $ \infty $ . The argument $ \frac{\omega_R}{\omega} $ inside the $ Li_2 $ function is greater than one. I do not much understand the $ Li_2 $ function, but I guess any number greater than one is outside its valid domain.

So, the major problem is the third part of the integral seems to not converge. But I'm sure it's result is a finite value. The full combined integral must result a finite value. Otherwise the transform couldn't exist. I tried a numerical integration of this integral with an Excel sheet, with very high frequencies, and it's seems to converge..., very slowly, the integral stabilizes with increasing frequency.

I got stuck here. What should I do?, how can I continue?, did I mess somewhere?

$\endgroup$
7
  • $\begingroup$ Why not interpolate the frequency-domain data into linear spaced one? $\endgroup$
    – ZR Han
    Feb 11, 2022 at 2:21
  • $\begingroup$ because linear spacing has very poor resolution in the low frequencies, with log spaced data you have good resolution at both extremes of the audio spectrum. in the audio frequency range, log frequency scale is used. $\endgroup$ Feb 11, 2022 at 2:26
  • $\begingroup$ I mean interpolate log scaled frequencies [2, 4, 8, 16, 32] into [2, 4, 6, 8, ... 32] and then you can apply FFT. Also you may consider nufft, which stands for nonuniform fast Fourier transform. $\endgroup$
    – ZR Han
    Feb 11, 2022 at 2:32
  • $\begingroup$ The software I want to copy performs this transform without using fft, in the user manual says it performs a properly integration from 0 to infinity. But I'll check that of nufft, sounds interesting, thanks. $\endgroup$ Feb 11, 2022 at 2:41
  • $\begingroup$ Seems a duplicate of electronics.stackexchange.com/questions/607025/…. Why do you say "FFT is not allowed"? FFT is the standard solution, especially for a system with zeros in the transfer function, which makes the integral improper. Also, you are using the definition of a continuous Hilbert transformer, not a discrete one. $\endgroup$
    – Hilmar
    Feb 11, 2022 at 14:16

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.