Формулы преобразования аналоговых фильтров в цифровые, Easyelectronics: Цифровые рекурсивные фильтры. Часть 1

import numpy as np

class ConvolutionDef_P2:
    def __init__(self):
        self.B = [1.0, 0.0, 0.0]
        self.A = [0.0, 0.0]

def HPF2_CalcScales(freqHz, sampleHz):
    if (freqHz > 0) and (sampleHz > 0):        
        try:
            betta  = np.tan(freqHz / sampleHz * np.pi)
            betta2 = betta * betta

            A0 = 1 + np.sqrt(2) * betta + betta2
            A1 = 2 - 2 * betta2
            A2 = 1 - np.sqrt(2) * betta + betta2
            B0 = betta2
            B1 = -2 * betta2
            B2 = betta2

            convol = ConvolutionDef_P2()
            convol.B[0] = B0 / A0
            convol.B[1] = B1 / A0
            convol.B[2] = B2 / A0
            convol.A[0] = -A1 / A0
            convol.A[1] = -A2 / A0

        except Exception as e:
            convol = Convolution()
            print("ERROR FilterCalc: {}".format(str(e)))

    print(vars(convol))
    return convol

def LPF2_CalcScales(freqHz, sampleHz):
    scales = HPF2_CalcScales(freqHz, sampleHz)
    scales.A[0] = -scales.A[0]
    scales.B[1] = -scales.B[1]
    return scales

def FilterApplyToSample(scales, sample):
    resp = []
    resp = np.zeros(len(sample))

    X = []
    X = np.zeros(2)
    Y = []
    Y = np.zeros(2)

    i = 0
    for x in sample:
        y =   x * scales.B[0] + X[0] * scales.B[1] + X[1] * scales.B[2] + \
           Y[0] * scales.A[0] + Y[1] * scales.A[1]
        X[1] = X[0]
        X[0] = x
        Y[1] = Y[0]
        Y[0] = y

        resp[i] = y
        i += 1

    return resp

Например, из суммы синусов с некоторым шагом, для подачи на вход фильтра:

def FillSampleByArr(sett, freqsArr):
    end_time = sett.samplePoints / sett.sampleRateHz
    time = np.arange(0, end_time, 1 / sett.sampleRateHz)

    sinewaves = np.zeros(len(time))
    n = 0
    for freq in freqsArr:
      sine_i = sett.sampleAmp * np.sin((2 * np.pi) * freq * time)
      sinewave += sine_i
      n += 1
    sinewaves = sinewaves / n
    sinewaves = sinewaves.astype(int)
    return sinewaves, time

def FillSampleByEvenFreq(sett, freqHzFrom, freqHzTo, freqHzStep):
    freqsArr = np.arange(freqHzFrom, freqHzTo, freqHzStep)
    return FillSampleByArr(sett, freqsArr)

import matplotlib.pyplot as plt

sett = SampleSettings(
    samplePoints = 2000,
    sampleRateHz = 50000,
    sampleAmp    = 0x01FF
)

fltCutFreqHz = 5000
    
sample, time = FillSampleByEvenFreq(sett, freqFromHz = 400, freqToHz = sett.sampleRateHz / 4, freqStepHz = 200)    

# Filter Select:
filterScales = HPF2_CalcScales(freqHz = fltCutFreqHz, sampleHz = sett.sampleRateHz)
#filterScales = LPF2_CalcScales(freqHz = fltCutFreqHz, sampleHz = sett.sampleRateHz) 

response = FilterApplyToSample(filterScales, sample)

# FFT
samplePoints = len(sample)
freqPoints = int(samplePoints/(2 * freqAxisDiv))        

# Frequency domain representation
fftSample = np.fft.fft(sample) / samplePoints  # Normalize amplitude
fftSample = fftSample[range(freqPoints)]       # Exclude sampling frequency
fftSample = abs(fftSample)

fftResponse = np.fft.fft(response) / samplePoints
fftResponse = fftResponse_dbg[range(freqPoints)]
fftResponse = abs(fftResponse)

# Show:
plt.figure()

plt.subplot(221)
plt.plot(time, sample)
plt.subplot(222)
plt.plot(time, response)

values      = np.arange(freqPoints)
timePeriod  = samplePoints / sett.sampleRateHz
axFreqs = values / timePeriod

plt.subplot(223)
plt.plot(axFreqs, fftSample)
plt.subplot(224)
plt.plot(axFreqs, fftResponse)

plt.show()