scipy.signal.minimumphase 偏离预期幅度

Question

我有一个问题，scipy.signal.minimumphase 没有提供与我输入相同的幅度响应，它偏离了很多。

长话短说。我有一个 material，其中所述 material 的吸收是在倍频程（6 个离散值）中测量的。我需要符合这些值的脉冲响应。我想通过使用最小相位 FIR 滤波器来近似脉冲响应，因为我没有任何相位信息。

我的程序如下

使用样条插值方案从离散值创建所需的幅值响应。
镜像内插幅度响应
IFFT 幅度响应。
将脉冲响应移动其长度的一半，创建线性相位并使其具有因果关系
使用scipy.signal.minimuphase创建最小相位脉冲

样条插值不是这里的问题，问题是最小相位的幅度偏离线性相位很多，我似乎找不到错误或原因？？运行下面的代码会产生这些响应：

这是代码的最小工作部分，它绘制了内插幅度响应、线性相位幅度和最小相位幅度。

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

我在做同样的事情，只是我使用scipy.signal.firwin2直接根据一些倍频程增益（你的步骤1~4）设计线性相位滤波器，然后使其成为最小相位。

如果您查看 documentation 中的参数详细信息，您会看到：

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

并且您的情节反映了这种平方根关系。简单的解决方法是在将 "desired magnitude response" 馈送到 scipy.signal.minimuphase.

之前预先平方

scipy.signal.minimumphase 偏离预期幅度

scipy.signal.minimumphase deviates from desired magnitude

filtering

signal-processing

scipy

digital-filter

python-3.x