Why extracted minimum phase component have inverted phase?

Question

I'm using Octave code to extract minimum phase, but this minimum phase has became inverted phase acording to original signal, What is goin wrong?

Resault here I have shifted main signal about 48 sample to matches it's minimumphase component(orange) and obviously these are out of phase:

Note:

The coefficients length is 2000 and the algorithm worked well for small length coefficients (15 weight tested).

The code:

function  [min,num_ap,den_ap,ap_zeroes,ap_poles] = minalldecompV2(b)
% impulse response must be given in vector b

tol=0.002;
  k = abs(b(1));                  % must be non-zero, gan factor
  zz = roots(b);                  % find zeros

  I = find( abs(zz) > 1+tol );
  zzr = zz(I);                    % zeros that need to be reflected

  scale = prod( abs( zzr ) );     % scaling factor for correct magnitude
  num_ap=poly(leja(zzr))/scale;
  ap_zeroes=zzr;
  zzr = 1 ./ conj(zzr);           % reflect zeros
  den_ap=poly(leja(zzr));
  ap_poles=zzr;
  min_z(I) = zzr;                         % replace zeros outside circle by reflected zeros
  bm = k * scale * poly( leja(min_z) )'; %leja    % find polynomial coefficients and apply scaling

if isreal( b )% remove (small) imaginary part if original filter is real-valued
  bm = real(bm); endif

min=hm;

I've found that the problem is general: Seems it act with absoluted version of signal: As you can see we have same minimumphase signal despite original signals have different polarity.

Answer 1

I believe the issue is with taking the absolute value of the first sample in the impulse response as the "gan" factor. If the first sample is negative, this will invert the result as implemented. Whether the resulting minimum phase response is inverted or not will come out of the process of simply reflecting zeros outside the unit circle to be inside the unit circle and computing the new polynomial for that result.

To prove this I will first review of "Min Phase" as well as the other common *-Phase systems for all causal stable discrete-time systems. The poles must be inside the unit circle for all, and it is the location of the zeros that distinguishes each:

Minimum phase system: All zeros inside the unit circle.

Maximum phase system: All zeros outside the unit circle.

Mixed phase system: Zeros both inside and outside the unit circle.

All pass system: All zeros at at the reciprocal ($1/p^*$) location of all poles.

Linear phase system All poles are at the origin and has reciprocal zeros inside and outside the unit circle.

Any arbitrary causal stable system that isn't already minimum phase can be decomposed into a minimum phase system cascaded with an all-pass system. This is depicted in the pole-zero diagrams below for an example mixed-phase system:

This was done by reflecting (conjugate inversion) the zeros to be inside the unit circle. What we see in the graphic is how the new poles for the all-pass system are exactly where the zeros are located in the minimum phase system. The pole zero cancellation at these locations when the systems are cascaded results in the original mixed-phase system.

The relationship given above mathematically is:

$$H(z)_{mix} = K\frac{(z-q_1)(z-q_1^*)(z-q_2)(z-q_2^*)}{(z-p_1)(z-p_1^*)z^2} $$

With $q_2$ and it's complex conjugate $q_2^*$ representing the two zeros outside the unit circle. $K$ is an arbitrary gain constant for the transfer function.

$$H(z)_{mix} = K_1K_2H(z)_{min}H(z)_{all} $$

With:

$$H(z)_{min} = K_1\frac{(z-q_1)(z-q_1^*)(z-1/q_2^*)(z-1/q_2)}{(z-p_1)(z-p_1^*)z^2} $$

$$H(z)_{all} = K_2\frac{(z-q_2)(z-q_2^*)}{(z-p_1)(z-p_1^*)} $$

The procedure to extract the minimum phase system given an arbitrary transfer function is simply:

• Find the roots
• Mirror the roots outside the unit circle using complex conjugate inversion
• Find the polynomial with the new roots, adjusting the gain using a location where the magnitude has a sufficient response. For example with a low pass adjust the gains so the frequency response at DC matches: scale $H_{min}$ by $H_{mix}(z=0)/H_{min}(z=0)$, for a high pass adjust the gains so the frequency response at Nyquist matches: scale $H_{min}$ by $H_{mix}(z=-1)/H_{min}(z=-1)$.

For very large polynomials, Leja ordering is done for numerical stability as MattL has referenced at this post here. Note that there is no test for the sign of the original transfer function. The coefficients for the minimum phase system may be inverted or not depending on the specific polynomials used, but this will come out of the process above directly. The coefficients will certainly be different unless the original system was already minimum phase, including sign inversions when necessary, but that does NOT mean the signal that passes through the system will be inverted. This depends on the frequency content of the signal and the resulting phase response at that frequency. I demonstrate this with the examples below.

Examples

Consider the simplest case of a linear phase filter with coefficients [-1, 2, -1] to test the procedure. This has roots complementary roots at $z=3.732$ and $z=0.2678$. The solution should simply reflect the root that is outside the unit circle so that both roots are at $z=0.2678$ resulting in coeff [1, -.535898, 0.0718].

Here is another example that is not linear phase along with showing the resulting frequency response where we see the "minimum phase" result:

The transfer function is given as:

$$H(z) = -1 + 4z^{-1} + z^{-2}$$

(The "coefficients" in this case are [-1, 4, 1])

The roots are [ 4.23606798 -0.23606798].

We reflect the first root to be inside the unit circle resulting in [ 0.23606798 -0.23606798]

Multiplying that out we get: [ 1, 0, -0.05572809]

We can use the DC response or $z=0$ in this case to match the scale, which results in a simple ratio of the sum of the coefficients for each to get [ 4.23606798, 0, -0.23606798]

So the minimum phase system with the exact same magnitude response but different phase would be as follows:

$$H(z)_{min} = 4.23606798 - 0.23606798z^{-2}$$

Below is a comparison of the frequency response for each where as expected we get the exact same magnitude response but the phase deviation over frequency has a minimum excursion (thus minimum phase), and the smallest derivative with respect to frequency (so minimum delay):

Note that both frequency responses of phase in this case both started at 0 degrees, yet the polynomial describing the transfer function did change sign (we went from [-1, 4, 1] to [4.236, 0, -0.236]). We didn't change this sign specifically but came out in the process of reflecting the zero and multiplying it out to get the result. If we were to change the sign based on the -1 coefficient, the result above would have 180 degrees added to it which would not be correct (although still a minimum phase system, but not the minimum phase equivalent).

What we see from the above plot is a signal that is close to DC will have the same phase at the input and output for both cases (original or minimum phase), but a signal that is close to Nyquist (half the sampling frequency) will be inverted between the input and output for the original system, but be in-phase for the minimum phase system.