scipy sosfilt 数组形状

Question

我正在尝试以两种方式绘制滤波器响应，一种是针对时间绘制，另一种是针对频率绘制。我现在的目标是绘制级联滤波器、filter1 和 filter2。 与频率的关系是我陷入困境的地方。 sosfreqz 期望 sos 数组具有 (n_sections, 6) 形状。 有好心人帮忙吗？我不知道如何解决这个问题，还没有在网上找到解决方案。 我让它以 (b, a) 格式工作，但 sos 显然在数字上更准确。

import numpy as np
from scipy.signal import butter, lfilter, freqs, sosfilt, sosfiltfilt, convolve, sosfreqz
import matplotlib.pyplot as plt

def butter_lowpass_filter(signal, freq_cutoff, sampling_freq, filter_order):
    sos = butter(filter_order, freq_cutoff, fs=sampling_freq , btype='low', analog=False, output='sos')
    return sosfilt(sos, signal )

#***********************************************************************
#*   Inputs - values rads/s converted to Hz at plot
#***********************************************************************
frequency = 10                            # Frequency, Hz
amplitude = 1                               # Amplitude
theta = 0
duration = 1                                # Duration of signal, seconds
sample_rate = 44100                         # Samples per second
no_of_samples = int(duration * sample_rate) # total number of samples

#***********************************************************************
#*   generate X-axis (time)
t = np.linspace(0, duration, no_of_samples, endpoint=False)

#***********************************************************************
#* generate white noise - Y-axis
noise = np.random.normal(0, 1, no_of_samples)
noise = amplitude * (noise / np.max(np.abs(noise)))  # Normalize the white noise
#noise = (noise * 2**15).astype(np.int16)           # Convert the white noise to a 16-bit format

#***********************************************************************
#*   generate sin wave - Y-axis
sinewave = amplitude * np.sin(2 * np.pi * frequency * t + theta) + np.sin(2*np.pi*20*t)    # rad/s to Hz

signal = sinewave       # use this to change between noise & sinewave

#***********************************************************************
#*   Subplot 1 - Freq vs Time plot, Plot the original signal.
plt.figure(1)
plt.subplot(2,1,1)
plt.plot(t, signal, label='input signal')
plt.xlabel('Time [sec]')
plt.grid()
plt.legend(loc='upper right')

#***********************************************************************
#*   Filter - stage 1
order = 2       # filter order number (20dB per decade)
cutoff = 27.5   # desired cutoff frequency of the filter, Hz

y1 = butter_lowpass_filter(signal, cutoff, sample_rate, order)

#***********************************************************************
#*   Filter - stage 2
order = 2       # filter order number (20dB per decade)
cutoff = 121   # desired cutoff frequency of the filter, Hz

y2 = butter_lowpass_filter(signal, cutoff, sample_rate, order)

#***********************************************************************
#*   combine filters
y = convolve(y1, y2, mode='same')   # combine arrays, but maintain original size 'same'

#***********************************************************************
#*   Subplot 2 - Freq vs Time plot, Plot filtered signal.
plt.subplot(2,1,2)
plt.plot(t, y, label='filtered_signal')
plt.xlabel('Time [sec]')
plt.grid()
plt.legend(loc='upper right')

y = y1          # use this to change between y, y1, y2 

print(len(y), type(y), np.shape(y), y)

w, h = sosfreqz(y, worN=44100)
plt.figure(2)
#plt.subplot(2,1,3)
#plt.semilogx(w, 20 * np.log10(abs(h)))              # Radians
plt.semilogx(w, 20 * np.log10(abs(h)) / (2* np.pi))  # Hz
plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [Hz]')
#plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(cutoff, color='green') # cutoff frequency

plt.show()

Answer 1

您可以使用

sosfiltfilt()

并将二阶部分组合到

sosfreqz

并绘制频率响应：

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


def butter_lowpass_filter(signal, cutoff_freq, sampling_freq, filter_order):
    sos = butter(filter_order, cutoff_freq, fs=sampling_freq, btype='low', analog=False, output='sos')
    return sosfiltfilt(sos, signal)


def _gen_signal(amp, freq, duration, rate):
    t = np.linspace(0, duration, duration * rate, endpoint=False)
    signal = amp * np.sin(2 * np.pi * freq * t) + np.sin(2 * np.pi * 20 * t)
    return t, signal


def _plot_time(t, signal, filtered_signal):
    plt.figure(1)
    plt.subplot(2, 1, 1)
    plt.plot(t, signal, label='input signal')
    plt.xlabel('Time [sec]')
    plt.grid()
    plt.legend(loc='upper right')

    plt.subplot(2, 1, 2)
    plt.plot(t, filtered_signal, label='filtered_signal')
    plt.xlabel('Time [sec]')
    plt.grid()
    plt.legend(loc='upper right')


def _plot_frequency(sos, cutoff_freq, sample_rate):
    w, h = sosfreqz(sos, worN=sample_rate)
    plt.figure(2)
    plt.semilogx(w, 20 * np.log10(abs(h)) / (2 * np.pi))
    plt.title('Butterworth filter frequency response')
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude [dB]')
    plt.margins(0, 0.1)
    plt.grid(which='both', axis='both')
    plt.axvline(cutoff_freq, color='green')
    plt.show()


def run(*args):
    amplitude, frequency, theta, duration, sample_rate = args
    t, signal = _gen_signal(amplitude, frequency, duration, sample_rate)

    order, cutoff = 2, 27.5
    sos1 = butter(order, cutoff, fs=sample_rate, btype='low', analog=False, output='sos')
    y1 = butter_lowpass_filter(signal, cutoff, sample_rate, order)

    order, cutoff = 2, 121
    sos2 = butter(order, cutoff, fs=sample_rate, btype='low', analog=False, output='sos')
    y2 = butter_lowpass_filter(signal, cutoff, sample_rate, order)

    stack = np.vstack((sos1, sos2))
    sosff = sosfiltfilt(stack, signal)

    _plot_time(t, signal, sosff)
    _plot_frequency(stack, cutoff, sample_rate)


if __name__ == "__main__":
    amp, freq = 1, 10
    theta, duration, sample_rate = 0, 1, 44100
    run(amp, freq, theta, duration, sample_rate)

打印