数字滤波器(4)—IIR/FIR系统对连续采集数据的滤波处理和模拟仿真

# IIR 系统分段处理：过去状态的处理
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from numpy.fft import fft, fftfreq

# 创建模拟信号
# 采样频率
Fs = 1000
nyq = 0.5 * Fs

# 设定阻带频率
low = 45.0
high = 55.0

# 按照通带和阻带频率的顺序来设计滤波器
lowcut = low / nyq
highcut = high / nyq

# 使用iirdesign设计滤波器
b, a = signal.iirdesign([lowcut, highcut], [lowcut - 0.05, highcut + 0.05], gpass=1, gstop=60, ftype='butter', output='ba')

T = 1.0 / Fs # 采样时间间隔
t = np.arange(0, 1, T) # 创建时间轴
x = np.sin(2 * np.pi * 50 * t) + 1 * np.sin(2 * np.pi * 220 * t) # 50Hz和120Hz信号叠加
# 创建一个与x长度相同的用于保存滤波后结果的数组
y = np.zeros_like(x)

# 分段进行滤波
win_size = 100 # 窗口大小
zi = signal.lfilter_zi(b, a) * x[0] # 计算初始状态
for i in range(win_size, len(x)+1, win_size):
    y[i-win_size:i], zi = signal.lfilter(b, a, x[i-win_size:i], zi=zi)

# 计算输入信号和输出信号的频谱
frequencies = fftfreq(len(t), 1/Fs)
x_spectrum = np.abs(fft(x))
y_spectrum = np.abs(fft(y))

# 绘制幅频响应图
w, h = signal.freqz(b, a)
plt.figure()
plt.plot(w/(2*np.pi), abs(h))
plt.title('Frequency response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Gain')
plt.grid()

# 绘制频谱图
plt.figure(figsize=(9, 6))
plt.subplot(221)
plt.title('Spectrum Before filtering')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')
plt.grid()
plt.plot(frequencies[:Fs//2], x_spectrum[:Fs//2], label='Before filtering')
plt.subplot(222)
plt.plot(frequencies[:Fs//2], y_spectrum[:Fs//2], label='After filtering')
plt.title('Spectrum After filtering')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')
plt.grid()

# 绘制模拟信号和滤波后的信号
plt.subplot(223)
plt.plot(t, x, label='Before filtering')
plt.title('Time domain Before filtering')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.grid()
plt.subplot(224)
plt.plot(t, y, label='After filtering')
plt.title('Time domain After filtering')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.grid()
plt.tight_layout()
plt.legend()

plt.show()

# FIR系统时域卷积：重叠保留法
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fft import fft,fftfreq

# 设置采样率和数据长度
fs = 6000

# 创建带通滤波器
# Create a bandpass filter
f1 = 400
f2 = 600
filter_len = 100         # FIR滤波器长度
b = signal.firwin(filter_len, [f1, f2], pass_zero=False, fs=fs, window='hamming')  #boxcar 

# 设置数据长度：模拟每次采集同样长度的数据
# segment_len，是按照np.concatenate函数在Valid模式下卷积后有效长度来计算的
seg_filter_len = 512                            # filter output length of each segment data
segment_len = seg_filter_len - filter_len + 1   # 分段数据目标长度 seg_filter_len = segment_len + filter_len - 1
target_length = segment_len * 50                # 总数据长度

# 而新的时间序列的上限
TimeSpace = target_length / fs

# 生成的时间序列为L的整数倍，模拟每次采样的数据的长度
t = np.arange(0, TimeSpace, 1/fs)

# 产生一个含有1.4KHz，2.4KHz和500Hz信号的模拟信号，其中1.4，2.4KHz的信号将被过滤
x = np.sin(2 * np.pi * 500 * t) + 0.5 * np.sin(2 * np.pi * 1400 * t) + 0.5 * np.sin(2 * np.pi * 2400 * t)

# Filter and segment with overlap-save method
# each segment length should be longer than filter's length
y_result = []
segment = None
for i in range(0, target_length, segment_len):
    if i==0:
        # 第一段数据滤波前，前面填充(M-1)个0，再加开始的L个数据参与滤波卷积
        segment = np.concatenate((np.zeros(filter_len-1), x[0:segment_len])) 
    else:
        # 对后续的每个段，包括从上段的末尾取(M-1)个点，再新取L个点，这里M指滤波器长度,L指拟当前新取的数据长度
        segment = np.concatenate((x[i-filter_len+1:i], x[i:i+segment_len]))
    
    # 每个段使用“valid”卷积模式，所以每次卷积后的有效数据长度是:L = L + 2(M-1) - 2(M-1) = segment_len
    conv = np.convolve(segment, b, mode='valid')
    # 与结果进行合并连接
    y_result = np.concatenate((y_result, conv))

# Frequency response of filter
w, h = signal.freqz(b, 1, fs=fs)
plt.figure()
plt.plot(w, abs(h))
plt.title('Frequency Response')     # 频率响应
plt.xlabel('Frequency [Hz]')        # 频率
plt.ylabel('Amplitude')             # 幅度
plt.grid()

# Plot the original signal and its FFT
n = len(x)
freq = fftfreq(n, 1/fs)
y = fft(x)

plt.figure(figsize=(9,6))
plt.subplot(221)
plt.plot(t[:500], x[:500])
plt.title('Original Signal')                # 原始信号
plt.grid()
plt.subplot(222)
plt.plot(freq[:n//2], np.abs(y[:n//2]*2/n)) # 频谱规范化输出
plt.title('FFT of Original Signal')         # 原始信号的FFT
plt.grid()

# Plot the filtered signal and its FFT
n = len(y_result)
freq = fftfreq(n, 1/fs)
y = fft(y_result)

plt.subplot(223)
plt.plot(t[:500], y_result[:500])
plt.title('Filtered Signal')                # 滤波后的信号
plt.grid()
plt.subplot(224)
plt.plot(freq[:n//2], np.abs(y[:n//2]*2/n)) # 频谱规范化输出
plt.title('FFT of Filtered Signal')         # 滤波后信号的FFT
plt.grid()

plt.tight_layout()
plt.show()