智能诊断综合实验

实验三

1.绘制频谱，分析干扰噪声

原始信号代码：

import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

dataFile_100 = r"..\A00100.mat"
data_100 = scio.loadmat(dataFile_100)
for val in data_100.keys():
    Y_100 = data_100[val]

Fs = 360
N = 2048
frq = np.arange(N)  # 频率数2048个数
half_x = frq[range(int(N / 2))]  # 取一半区间
print("心电信号数据:", Y_100)
print("数据长度:", len(Y_100))

# 只取3000组数据
Y_100 = Y_100[0][:3000]

# 原始数据
plt.figure(figsize=(15, 12))
plt.subplot(111)
plt.title('原始心电信号数据')
plt.plot(range(0, len(Y_100)), Y_100)
plt.xlabel("N")
plt.ylabel("心电信号")
plt.show()

2.采用巴特沃斯滤波器

源码如下：

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

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

dataFile_100 = r"..\A00100.mat"
data_100 = scio.loadmat(dataFile_100)
for val in data_100.keys():
    Y_100 = data_100[val]

Fs = 360

# 原始数据
plt.figure(figsize=(15, 12))
plt.subplot(311)
plt.title('原始心电信号数据')
plt.plot(range(0, len(Y_100[0])), Y_100[0])
plt.xlabel("N")
plt.ylabel("心电信号")

# 低通滤波器设计和滤波
b, a = signal.butter(4, 100, 'low', analog=True)  # 4阶低通临界频率为100Hz
w, h = signal.freqs(b, a)  # h为频率响应,w为频率
result_lowpass = signal.lfilter(b, a, Y_100[0])

plt.subplot(312)
plt.title('低通滤波后的心电信号数据')
plt.plot(range(0, len(result_lowpass)), result_lowpass)
plt.xlabel("N")
plt.ylabel("心电信号（低通滤波后）")

plt.show()

结果：

巴特沃斯带通：

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

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

dataFile_100 = r"..\A00100.mat"
data_100 = scio.loadmat(dataFile_100)
for val in data_100.keys():
    Y_100 = data_100[val]

Fs = 360

# 原始数据
plt.figure(figsize=(15, 12))
plt.subplot(311)
plt.title('原始心电信号数据')
plt.plot(range(0, len(Y_100[0])), Y_100[0])
plt.xlabel("N")
plt.ylabel("心电信号")

# 带通滤波器设计和滤波
N, Wn = signal.buttord([20, 50], [14, 60], 3, 40, True)
b, a = signal.butter(N, Wn, 'band', True)
w, h = signal.freqs(b, a, np.logspace(1, 2, 500))
result_bandpass = signal.lfilter(b, a, Y_100[0])

plt.subplot(312)
plt.title('带通滤波后的心电信号数据')
plt.plot(range(0, len(result_bandpass)), result_bandpass)
plt.xlabel("N")
plt.ylabel("心电信号（带通滤波后）")

plt.show()

运行结果：

3.当阶数不同时运行结果也不同，N=8

N=6

实验四

1.设计的滤波器是巴特沃斯低通滤波器

代码如下：

# 设计低通滤波器去除100Hz以上心电干扰
cutoff_freq = 100  # 截止频率设为100Hz
b, a = signal.butter(4, cutoff_freq / (Fs / 2), 'low')  # 4阶低通滤波器设计
filtered_signal = signal.filtfilt(b, a, Y_100_first_3000)

2.中值滤波

源码：

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

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

dataFile_100 = r"..\A00100.mat"
data_100 = scio.loadmat(dataFile_100)
for val in data_100.keys():
    Y_100 = data_100[val]

Y_100_first_100 = Y_100[0][:3400]  # 仅取前100组数据

Fs = 360
N = 3400  # 仅处理前100组数据
frq = np.arange(N)  # 频率数100个数
half_x = frq[range(int(N / 2))]  # 取一半区间
print("心电信号数据:", Y_100_first_100)
print("数据长度:", len(Y_100_first_100))

# 原始数据（前100组）
plt.figure(figsize=(15, 12))
plt.subplot(311)
plt.title('原始心电信号数据（前100组）')
plt.plot(range(0, len(Y_100_first_100)), Y_100_first_100)
plt.xlabel("N")
plt.ylabel("心电信号")

# IIR零相移数字滤波器设计和滤波
Wp = 1.4 * 2 / Fs  # 通带截止频率
Ws = 0.6 * 2 / Fs  # 阻带截止频率
devel = 0.005  # 通带纹波
Rp = 20 * np.log10((1 + devel) / (1 - devel))  # 通带纹波系数
Rs = 20  # 阻带衰减

N, Wn = signal.ellipord(Wp, Ws, Rp, Rs, True)  # 求椭圆滤波器的阶次
b, a = signal.ellip(N, Rp, Rs, Wn, 'high', True)  # 求椭圆滤波器的系数
w, h = signal.freqs(b, a, 512)
result = signal.lfilter(b, a, Y_100_first_100)

plt.subplot(312)
plt.title('线性滤波后信号（前3000组）')
plt.plot(range(0, len(result)), result)
plt.xlabel("N")
plt.ylabel("心电信号（滤波后）")

# 绘制频谱图
N_fft = 512
plt.subplot(313)
plt.plot(np.abs(np.fft.fft(Y_100_first_100)) * 2 / N_fft)
plt.title('原始心电信号频谱（前3000组）')
plt.xlabel("(Hz)")
plt.ylabel("幅度")

plt.tight_layout()
plt.show()

运行结果：

3.巴特沃斯滤波器完整代码

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

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

dataFile_100 = r"..\A00100.mat"
data_100 = scio.loadmat(dataFile_100)
for val in data_100.keys():
    Y_100 = data_100[val]

Y_100_first_3000 = Y_100[0][:3000]  # 仅取前3000组数据

Fs = 360
N = 3000  # 仅处理前3000组数据
frq = np.arange(N)  # 频率数3000个数
half_x = frq[range(int(N / 2))]  # 取一半区间
print("心电信号数据:", Y_100_first_3000)
print("数据长度:", len(Y_100_first_3000))

# 原始数据（前3000组）
plt.figure(figsize=(15, 12))
plt.subplot(311)
plt.title('原始心电信号数据（前3000组）')
plt.plot(range(0, len(Y_100_first_3000)), Y_100_first_3000)
plt.xlabel("N")
plt.ylabel("心电信号")

# 设计低通滤波器去除100Hz以上心电干扰
cutoff_freq = 100  # 截止频率设为100Hz
b, a = signal.butter(4, cutoff_freq / (Fs / 2), 'low')  # 4阶低通滤波器设计
filtered_signal = signal.filtfilt(b, a, Y_100_first_3000)
plt.subplot(312)
plt.title('去除100Hz以上心电干扰后的心电信号数据（前3000组）')
plt.plot(range(0, len(filtered_signal)), filtered_signal)
plt.xlabel("N")
plt.ylabel("心电信号")

plt.show()

运行结果：