详细阐述ECG 中IIR滤波与FIR滤波
时间: 2023-03-14 17:35:30 浏览: 281
ECG中的IIR滤波和FIR滤波是两种不同的数字滤波技术,其中IIR滤波通过利用模拟滤波器中的模拟滤波器元件,以及有限低通,高通,带通或椭圆滤波器来实现,而FIR滤波则是一种线性滤波技术,它通过利用延迟线路和系数来构建理想的频率响应来实现。
相关问题
给定一组干净心电信号数据,数据文件存于C盘Ecg.txt。采样频率f= 500Hz。1、编写程序读出心电信号,并在屏幕上打印出其波形。 2、产生模拟高斯白噪声信号,与干净心电混合,设计一个IIR低通滤波器和一个FIR低通滤波器分别滤除心电信号中的白噪声干扰,调整白噪声信噪比大小,对滤波前后的心电信号的频谱进行分析比较。其中数字低通滤波器指标要求,通带截止频率Wp=0.1 Π,阻带截止频率Ws=0.16Π,阻带衰减不小于15 dB,通带衰减不大于1 dB。 要求:编写一个IIR低通滤波器和一个FIR低通滤波器仿真程序,在屏幕上打印出数字滤波器的频率区间[0,Π]上的幅频响应特性由线(H(ejw));计算其对含噪心电信号的低通滤波响应序列,并在屏幕上打印出干净心电信号波形,含工频干扰的心电信号波形以及IIR低通滤波和FIR低通后的信号波形,并进行比较;同时对滤波前后的心电信号的频谱进行分析比较,并在屏幕上打印出滤波前后的心电信号的频谱,观察其变化。
很抱歉,前面回答的代码有些问题,这里提供更正的Python代码供参考:
```python
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 500 # 采样频率
# 读取心电信号并打印波形
with open('C:/Ecg.txt', 'r') as f:
data = f.readlines()
data = np.array([float(d) for d in data]) # 转换为numpy数组
t = np.linspace(0, len(data)/fs, len(data)) # 时间序列
plt.plot(t, data)
plt.xlabel('Time (s)')
plt.ylabel('Voltage (mV)')
plt.title('Clean ECG Signal')
plt.show()
# 产生高斯白噪声信号并混合
def generate_gaussian_noise(sigma, length):
return np.array([np.random.normal(0, sigma) for i in range(length)])
noise = generate_gaussian_noise(sigma=0.1, length=len(data))
mixed_signal = data + noise
# IIR低通滤波器设计
Wp = 0.1*np.pi
Ws = 0.16*np.pi
Rp = 1
As = 15
b, a = signal.iirdesign(Wp, Ws, Rp, As, analog=False, ftype='cheby2', fs=fs)
# FIR低通滤波器设计
N = 101
Wn = 0.1*np.pi
b_fir = signal.firwin(N, Wn, window='hamming', pass_zero=True, fs=fs)
# 计算滤波器频率响应并打印
w, h_iir = signal.freqz(b, a)
w, h_fir = signal.freqz(b_fir, 1)
plt.plot(w/np.pi*fs/2, np.abs(h_iir), label='IIR')
plt.plot(w/np.pi*fs/2, np.abs(h_fir), label='FIR')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('Digital Lowpass Filter Frequency Response')
plt.legend()
plt.show()
# IIR低通滤波
filtered_signal_iir = signal.filtfilt(b, a, mixed_signal)
# FIR低通滤波
filtered_signal_fir = signal.filtfilt(b_fir, 1, mixed_signal)
# 打印心电信号波形和滤波后波形
plt.plot(t, data, label='Clean ECG')
plt.plot(t, mixed_signal, label='Noisy ECG')
plt.plot(t, filtered_signal_iir, label='IIR Filtered ECG')
plt.plot(t, filtered_signal_fir, label='FIR Filtered ECG')
plt.xlabel('Time (s)')
plt.ylabel('Voltage (mV)')
plt.title('ECG Signal with Gaussian White Noise')
plt.legend()
plt.show()
# 打印滤波前后心电信号频谱
f, Pxx_den_clean = signal.periodogram(data, fs)
f, Pxx_den_mixed = signal.periodogram(mixed_signal, fs)
f, Pxx_den_iir = signal.periodogram(filtered_signal_iir, fs)
f, Pxx_den_fir = signal.periodogram(filtered_signal_fir, fs)
plt.semilogy(f, Pxx_den_clean, label='Clean ECG')
plt.semilogy(f, Pxx_den_mixed, label='Noisy ECG')
plt.semilogy(f, Pxx_den_iir, label='IIR Filtered ECG')
plt.semilogy(f, Pxx_den_fir, label='FIR Filtered ECG')
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD (V^2/Hz)')
plt.title('ECG Signal PSD')
plt.legend()
plt.show()
```
这段代码可以实现对心电信号的读取、高斯白噪声的产生和混合、IIR和FIR数字低通滤波器的设计和滤波、滤波前后心电信号的波形和频谱比较等功能。需要注意的是,产生高斯白噪声时应使用`numpy.random.normal`函数,而不是`random.gauss`函数。另外,FIR滤波器的设计应使用`signal.firwin`函数,而不是`signal.fir_design`函数。
clear,clc; val=importdata('Ecg.txt'); signal=val(1,1:1800); fs=500; figure(1) subplot(4,2,1); plot(signal); title('干净的EGC信号'); xlabel('采样点'); ylabel('幅值(dB)'); grid on; signal1=awgn(signal,10,'measured'); subplot(4,2,2); plot(signal1); title('高斯噪声的EGC信号'); xlabel('采样点'); ylabel('幅值(dB)'); % 设计IIR低通滤波器 Wp = 0.1pi; % 通带截止频率 Ws = 0.16pi; % 阻带截止频率 Rp = 1; % 通带衰减 Rs = 15; % 阻带衰减 [n, Wn] = buttord(Wp, Ws, Rp, Rs, 's'); [b, a] = butter(n, Wn); % 绘制数字低通滤波器的幅频响应 [H, w] = freqz(b, a, 512); f = w/pi500; subplot(4,2,3); plot(w/pi,20log10(abs(H))); xlabel('频率'); ylabel('幅值(dB)'); title('IIR低通滤波器幅频响应'); iir_filtered_signal = filter(b, a, signal1); subplot(4,2,4); plot(iir_filtered_signal); xlabel('频率'); ylabel('幅值(dB)'); title('IIR低通滤波后的含高斯噪声的图像'); iir_signal = abs(fft(signal)); subplot(4,2,5); plot(20log10(abs(iir_signal))); xlabel('频率'); ylabel('幅值(dB)'); title('含高斯噪声的频谱'); iir_signal1 = abs(fft(signal1)); subplot(4,2,6); plot(20log10(abs(iir_signal1))); xlabel('频率'); ylabel('幅值(dB)'); title('IIR低通滤波后的含高斯噪声的频谱'); n = 80; % 滤波器阶数 wc = 0.1pi; % 通带截止频率 h = fir1(n, wc/(fs/2), kaiser(n+1, 6)); % 计算FIR低通滤波器系数 filtered_signal_fir = filter(h, 1, signal); % 应用FIR滤波器 subplot(4,2,7); plot(20log10(abs(h))); title('FIR低通滤波幅频响应'); xlabel('频率'); ylabel('幅值(dB)'); [Pxx_filtered_fir, f_filtered_fir] = periodogram(filtered_signal_fir, [], [], fs); subplot(4,2,8); plot(20*log10(abs(Pxx_filtered_fir))); title('FIR低通滤波后的含高斯噪声的频谱'); xlabel('频率'); ylabel('幅值(dB)');逐句注释这段代码
这段代码是用 MATLAB 实现了对 ECG 信号进行数字滤波的过程。下面是逐句的注释:
```
clear,clc; % 清空命令窗口和工作空间
val=importdata('Ecg.txt'); % 读入 ECG 数据
signal=val(1,1:1800); % 取 ECG 数据的前 1800 个采样点
fs=500; % 采样率为 500 Hz
figure(1) % 创建一个新图形窗口
subplot(4,2,1); % 将窗口分为 4 行 2 列,并在第 1 个子图中绘制
plot(signal); % 绘制干净的 ECG 信号
title('干净的EGC信号'); % 设置标题
xlabel('采样点'); % 设置 x 轴标签
ylabel('幅值(dB)'); % 设置 y 轴标签
grid on; % 显示网格线
signal1=awgn(signal,10,'measured'); % 生成一个高斯噪声为 10 dB 的含噪声 ECG 信号
subplot(4,2,2); % 在第 2 个子图中绘制含噪声的 ECG 信号
plot(signal1); % 绘制含噪声的 ECG 信号
title('高斯噪声的EGC信号'); % 设置标题
xlabel('采样点'); % 设置 x 轴标签
ylabel('幅值(dB)'); % 设置 y 轴标签
% 设计 IIR 低通滤波器
Wp = 0.1*pi; % 通带截止频率
Ws = 0.16*pi; % 阻带截止频率
Rp = 1; % 通带衰减
Rs = 15; % 阻带衰减
[n, Wn] = buttord(Wp, Ws, Rp, Rs, 's'); % 计算滤波器的阶数和截止频率
[b, a] = butter(n, Wn); % 计算 IIR 低通滤波器的系数
% 绘制数字低通滤波器的幅频响应
[H, w] = freqz(b, a, 512); % 计算数字滤波器的频率响应
f = w/pi*500; % 将频率转换为 Hz
subplot(4,2,3); % 在第 3 个子图中绘制幅频响应
plot(w/pi,20*log10(abs(H))); % 绘制幅频响应
xlabel('频率'); % 设置 x 轴标签
ylabel('幅值(dB)'); % 设置 y 轴标签
title('IIR低通滤波器幅频响应'); % 设置标题
iir_filtered_signal = filter(b, a, signal1); % 对含噪声的 ECG 信号进行 IIR 低通滤波处理
subplot(4,2,4); % 在第 4 个子图中绘制滤波后的信号
plot(iir_filtered_signal); % 绘制滤波后的信号
xlabel('频率'); % 设置 x 轴标签
ylabel('幅值(dB)'); % 设置 y 轴标签
title('IIR低通滤波后的含高斯噪声的图像'); % 设置标题
iir_signal = abs(fft(signal)); % 计算干净的 ECG 信号的频谱
subplot(4,2,5); % 在第 5 个子图中绘制频谱
plot(20*log10(abs(iir_signal))); % 绘制幅值谱
xlabel('频率'); % 设置 x 轴标签
ylabel('幅值(dB)'); % 设置 y 轴标签
title('含高斯噪声的频谱'); % 设置标题
iir_signal1 = abs(fft(signal1)); % 计算含噪声的 ECG 信号的频谱
subplot(4,2,6); % 在第 6 个子图中绘制频谱
plot(20*log10(abs(iir_signal1))); % 绘制幅值谱
xlabel('频率'); % 设置 x 轴标签
ylabel('幅值(dB)'); % 设置 y 轴标签
title('IIR低通滤波后的含高斯噪声的频谱'); % 设置标题
n = 80; % 滤波器阶数
wc = 0.1*pi; % 通带截止频率
h = fir1(n, wc/(fs/2), kaiser(n+1, 6)); % 计算 FIR 低通滤波器的系数
filtered_signal_fir = filter(h, 1, signal); % 对干净的 ECG 信号进行 FIR 低通滤波处理
subplot(4,2,7); % 在第 7 个子图中绘制幅频响应
plot(20*log10(abs(h))); % 绘制幅频响应
title('FIR低通滤波幅频响应'); % 设置标题
xlabel('频率'); % 设置 x 轴标签
ylabel('幅值(dB)'); % 设置 y 轴标签
[Pxx_filtered_fir, f_filtered_fir] = periodogram(filtered_signal_fir, [], [], fs); % 计算 FIR 低通滤波后的信号的频谱
subplot(4,2,8); % 在第 8 个子图中绘制频谱
plot(20*log10(abs(Pxx_filtered_fir))); % 绘制幅值谱
title('FIR低通滤波后的含高斯噪声的频谱'); % 设置标题
xlabel('频率'); % 设置 x 轴标签
ylabel('幅值(dB)'); % 设置 y 轴标签
```
注:原代码中的 `%` 应改为 `pi`。
阅读全文