scipy.signal.minimumphase deviates from desired magnitude

Question

I have a problem with the scipy.signal.minimumphase not providing the same magnitude response as i input, it deviates quite a lot.

Long story short. I have a material, where the absorption of said material is measured in octave bands(6 discrete values). I need an impulse response that fit those values. I want to approximate that impulse response by using minimum phase FIR filters since i don't have any phase information.

My procedure is the following

1. Use a spline interpolation scheme to create a desired magnitude response from the discrete values.
2. Mirror the interpolated magnitude response
3. IFFT the magnitude response.
4. Shift the impulse response by half its length, creating linear phase and making it causal
5. Use scipy.signal.minimuphase to create a minimum phase impulse

The spline interpolation is not the issue here, the problem is that the magnitude of the minimum phase deviates alot from the linear phase and i can't seem to find the error or reason why?? Running the code below yields you these responses:

Here is a minimum working part of the code, it plots the interpolated magnitude response, the linear phase magnitude and the minimum phase magnitude.

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

plt.close('all')

## Determining Gradient between points
def slope(x1, y1, x2, y2):
    m = (y2-y1)/(x2-x1)
    return m


## Parameters for designing the frequency response
fs = 48000
N = 4096
alpha = np.array([1,0.07,0.17,0.40,0.55,0.65,0.65,1])
beta = np.sqrt(1-alpha)
f = np.array([0,125,250,500,1000,2000,4000,(fs/2)])

### Spline interpolation
faxis = np.linspace(0,fs/2,np.int64(N))
gradient = np.zeros((np.size(f),1))
x = np.zeros([4,np.size(f)-1])   

for n in range(0,np.size(f)-1,1):
    gradient[n] = slope(f[n],beta[n],f[n+1],beta[n+1])

for n in range(0,(np.size(f))-1,1):
    a = np.array([[f[n]**3,f[n]**2,f[n]**1,1], 
                  [f[n+1]**3,f[n+1]**2,f[n+1]**1,1],
                  [3*(f[n]**2),f[n]*2,1,0],
                  [3*(f[n+1]**2),f[n+1]*2,1,0]])
    b = np.array([beta[n],beta[n+1],gradient[n],gradient[n+1]])    
    x[:,n] = np.linalg.solve(a, b)

## Using a,b,c,d coef. to make the polynomials and then crop them to fit
poly = np.zeros([(np.size(faxis)),np.size(beta)])
combined = np.array([])
for n in range(0,np.size(poly,1)-1,1):
    poly[:,n] = x[0,n]*faxis**3+x[1,n]*faxis**2+x[2,n]*faxis+x[3,n]        
    combined = np.append(combined,poly[(faxis >= f[n]) & (faxis <= f[n+1]),n])        

## Mirror the single sided frequency response
Full = np.concatenate((combined,np.flipud(combined)))

## ifft the mirrored resposne - No Phase info
impulse = np.real(np.fft.ifft(Full))

## Add linear phase to the impulse - i.e. shift the response and make it causal 
impulseshift = np.roll(impulse,np.int(len(impulse)/2)-1)

## Change to minimumphase
minphase = signal.minimum_phase(impulseshift[:-1],method='hilbert',n_fft=fs)   


## Plot the results
wLinear, hLinear = signal.freqz(impulseshift,worN=4096)  
w, h = signal.freqz(minphase,worN=4096)  

plt.figure
plt.semilogx(faxis,combined,label="No Phase")
plt.semilogx((fs * 0.5 / np.pi) * w, abs(h),label="Minimum Phase")
plt.semilogx((fs * 0.5 / np.pi) * wLinear, abs(hLinear),label="Linear Phase")
plt.grid()
plt.legend()
plt.xlabel("Frequency [Hz]")
plt.ylabel("Amplitude [.]")

Answer 1

I'm doing the same thing except that I use scipy.signal.firwin2 to directly design the linear phase filter based on some octave band gains (your step 1~4) and then make it minimum phase.

If you look at the parameter details in the documentation you'll see:

the resulting minimum phase filter will have a magnitude response that approximates the square root of the original filter’s magnitude response.

And your plot reflects this square root relationship. The simple fix would be to pre-square your "desired magnitude response" before feeding them to scipy.signal.minimuphase.