Scipy sosfreqz报错：需(n_sections,6)形状sos数组绘制级联滤波器频响

问题核心 #

需要绘制两级级联Butterworth低通滤波器的时域、频域响应，时域响应正常，但调用sosfreqz时触发形状不匹配错误——该函数要求输入的sos数组必须为(n_sections,6)格式。此前用(b,a)格式可正常绘制频响，但sos格式数值精度更高，需解决该问题。

错误原因 #

1. 调用sosfreqz时传入了滤波后的信号数组，而非滤波器的sos系数矩阵，这是核心错误。
2. 原代码用convolve(y1, y2, mode='same')合并结果，这不是真正的滤波器级联（级联应是前级输出作为后级输入），且未保存两级滤波器的sos参数，无法计算级联后的频响。

解决方案步骤 #

1. 修改滤波器函数，同时返回滤波后的信号和对应的sos系数，保留每一级的滤波器参数。
2. 正确实现级联滤波：将第一级的输出作为第二级滤波器的输入。
3. 合并两级的sos系数矩阵：级联滤波器的sos数组只需将两级的sos数组垂直拼接（每一级对应一个shape为(1,6)的section，两级拼接后为(2,6)，符合(n_sections,6)要求）。
4. 使用合并后的sos系数调用sosfreqz，并利用fs参数直接获取Hz单位的频率轴，无需手动转换。

修正后的完整代码 #

import numpy as np
from scipy.signal import butter, sosfilt, sosfreqz
import matplotlib.pyplot as plt

def butter_lowpass_filter(signal, freq_cutoff, sampling_freq, filter_order):
    # 同时返回sos系数和滤波后的信号
    sos = butter(filter_order, freq_cutoff, fs=sampling_freq, btype='low', analog=False, output='sos')
    filtered_signal = sosfilt(sos, signal)
    return sos, filtered_signal

# 输入参数
frequency = 10                # 信号频率，Hz
amplitude = 1                 # 信号幅值
duration = 1                  # 信号时长，秒
sample_rate = 44100           # 采样率，Hz
no_of_samples = int(duration * sample_rate)

# 生成时间轴
t = np.linspace(0, duration, no_of_samples, endpoint=False)

# 生成测试信号：10Hz + 20Hz正弦波
sinewave = amplitude * np.sin(2 * np.pi * frequency * t) + np.sin(2*np.pi*20*t)
signal = sinewave

# 绘制原始信号
plt.figure(1)
plt.subplot(2,1,1)
plt.plot(t, signal, label='输入信号')
plt.xlabel('时间 [秒]')
plt.grid()
plt.legend(loc='upper right')

# 第一级滤波器：2阶，截止27.5Hz
order1 = 2
cutoff1 = 27.5
sos1, y1 = butter_lowpass_filter(signal, cutoff1, sample_rate, order1)

# 第二级滤波器：2阶，截止121Hz（级联：将y1作为输入）
order2 = 2
cutoff2 = 121
sos2, y2 = butter_lowpass_filter(y1, cutoff2, sample_rate, order2)

# 合并两级sos系数，得到级联滤波器的sos矩阵
sos_total = np.vstack([sos1, sos2])

# 绘制级联滤波后的时域信号
plt.subplot(2,1,2)
plt.plot(t, y2, label='级联滤波后信号')
plt.xlabel('时间 [秒]')
plt.grid()
plt.legend(loc='upper right')

# 绘制频域响应
plt.figure(2)
# 使用合并后的sos_total调用sosfreqz，指定fs参数直接获取Hz频率轴
w, h = sosfreqz(sos_total, worN=44100, fs=sample_rate)
plt.semilogx(w, 20 * np.log10(abs(h)))
plt.title('Butterworth级联滤波器频域响应')
plt.xlabel('频率 [Hz]')
plt.ylabel('幅值 [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
# 标注两级截止频率
plt.axvline(cutoff1, color='green', label=f'截止频率1: {cutoff1}Hz')
plt.axvline(cutoff2, color='red', label=f'截止频率2: {cutoff2}Hz')
plt.legend()

plt.show()

关键修改说明 #

• 函数butter_lowpass_filter新增返回sos系数，保留滤波器参数。
• 级联滤波改为sosfilt(sos2, y1)，符合实际滤波流程。
• 用np.vstack([sos1, sos2])合并两级sos，得到符合要求的(2,6)形状数组。
• sosfreqz调用时指定fs=sample_rate，直接输出Hz单位的频率值，无需手动转换。

内容的提问来源于stack exchange，提问作者lambcutlet