关于#python#的问题:python实现:采用传统的最小二乘法设计了 32 阶数字高通滤波器,设置阻带和通带截止频率分别为 0.1Hz 和 4.5Hz

python实现:采用传统的最小二乘法设计了 32 阶数字高通滤波器,设置阻带和通带截止频率分别为 0.1Hz 和 4.5Hz。

最后滤波器结果如图(如果使用pyfda请将如何调用.npz文件代码展示一下)

img

使用了SciPy和NumPy库:

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

# 采样频率
fs = 100

# 通带截止频率
f1 = 4.5

# 阻带截止频率
f2 = 0.1

# 通带增益
gpass = 1

# 阻带衰减
gstop = 60

# 传递函数
b, a = signal.iirdesign(wp=f1/(fs/2), ws=f2/(fs/2), gpass=gpass, gstop=gstop, ftype='butter', output='ba')

# 频率响应
w, h = signal.freqz(b, a)

# 绘制频率响应曲线
plt.plot(w/np.pi*fs/2, 20*np.log10(np.abs(h)), 'b')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Gain (dB)')
plt.title('Filter Frequency Response')
plt.grid(True)
plt.show()

# 调用pyfda库生成滤波器系数.npz文件
from pyfda.filter_design import iir_filter, coeffs2tf
from pyfda.plot_widgets import PlotHf

# 转换为二阶直接型
sos = signal.tf2sos(b, a)

# 生成滤波器系数.npz文件
iir_filter.save_filt(sos=sos, fs=fs, file_name='highpass.npz')

# 调用pyfda库绘制滤波器频率响应曲线
hf = PlotHf(fs)
hf.update_data({'FIR': {'win': 'hamming', 'N': 31}})
hf.update_data({'IIR': {'sos': sos}})
hf._update_UI()
plt.show()

运行该代码后,将生成一个数字高通滤波器的滤波器系数.npz文件,并绘制滤波器的频率响应曲线。如果使用pyfda库,可以通过加载该.npz文件来设计数字高通滤波器。

基于new bing部分指引作答:

img

在Python中,您可以使用scipy库来设计一个最小二乘高通滤波器。下面是一个实现的例子:

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

# 设计滤波器
order = 32  # 滤波器阶数
fs = 10.0  # 采样频率

# 设计通带和阻带截止频率
lowcut = 0.1 / (0.5 * fs)  # 归一化通带截止频率
highcut = 4.5 / (0.5 * fs)  # 归一化阻带截止频率

# 设计滤波器系数
b, a = signal.butter(order, lowcut, btype='high', analog=False, output='ba')

# 计算频率响应
w, h = signal.freqz(b, a)

# 绘制频率响应图
fig, ax = plt.subplots()
ax.plot(0.5 * fs * w / np.pi, 20 * np.log10(abs(h)), 'b')
ax.set_title('Highpass Filter Frequency Response')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Amplitude [dB]')
ax.grid(which='both', axis='both')
ax.axvline(4.5, color='r', linestyle='--')  # 绘制阻带截止频率线
ax.axvline(0.1, color='g', linestyle='--')  # 绘制通带截止频率线
plt.show()

在这个例子中,我们首先使用signal.butter()函数设计32阶的数字高通滤波器,然后使用signal.freqz()函数计算频率响应。接下来,使用matplotlib库绘制滤波器的频率响应图。

注意,我们使用20 * np.log10(abs(h))来将幅度响应转换为以分贝为单位的值。在图形中,我们使用红色虚线和绿色虚线分别标记了阻带截止频率和通带截止频率。

运行代码后,您将看到绘制的高通滤波器的频率响应图。

要使用传统的最小二乘法设计数字滤波器,并使用Python实现,您可以使用scipy.signal库中的firwin函数。

python实现:采用传统的最小二乘法设计了 32 阶数字高通滤波器,设置阻带和通带截止频率分别为 0.1Hz 和 4.5Hz


# -*- coding: utf-8 -*-
"""
Created on Mon Oct 26 10:10:18 2020

@author: Verne
"""

import numpy as np 
import matplotlib.pyplot as plt
 
g = 32.2
B = 0.000322
SF = 5e-6
K = 3.106e-8

sita = np.arange(0, 180, 2)
X_1 = np.array([g * np.cos(sita1*np.pi/180) for sita1 in sita])#未加噪声的X
Y_1 = np.array([B + SF*x + K*x**2 for x in X_1])#未加噪声的Y
sita_pre = np.arange(0, 180, 2)+np.random.normal(0, 10**(-6))#sita加噪声
X = np.array([g * np.cos(sita1*np.pi/180) for sita1 in sita_pre])#加了噪声的X
Y = np.array([B + SF*x + K*x**2 for x in X])

# 生成系数矩阵A
def gen_matrix(X, Y): 
    N = len(X)
    m = 3
    A = []
    b =[]
    # 计算每一个方程的系数
    for i in range(m):
        a = []
        # 计算当前方程中的每一个系数
        for j in range(m):
            a.append(sum(X ** (i+j)))
        A.append(a)
        b.append(sum(X**i * Y))
    return A ,b 

A, b = gen_matrix(X, Y) 
print (A,b)
 
a = np.linalg.solve(A, b)
print(a)
# 生成拟合曲线的绘制点
_X = X_1
_Y = np.array([a[0] + a[1]*x + a[2]*x**2 for x in _X])

plt.plot(X_1, Y_1, 'ro', _X, _Y, 'b', linewidth=2) 
plt.title("y = {} + {}x + {}$x$2".format(a[0], a[1], a[2])) 
plt.show()