数字滤波器（5）—FIR连续采样分段卷积时域重叠相加法

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

# 创建带通滤波器
f1 = 40
f2 = 60
filter_len = 80     # 滤波器长度
fs = 1600           # 采样频率维持不变

b = signal.firwin(filter_len, [f1, f2], pass_zero=False, fs=fs)

# 设置数据长度
seg_filter_len = 256                            # 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                # 总数据长度

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

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

# 产生一个含有300Hz,450Hz和50Hz信号的模拟信号
x = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 300 * t) + 0.5 * np.sin(2 * np.pi * 450 * t)

segments = []
for i in range(0, len(x), segment_len):
    segments.append(x[i:i+segment_len])

# Filtering & Overlap-Add processing
# Total outputput buffer, len = target_length + filter_len - 1
filtered_signals = np.zeros(target_length + filter_len - 1)             
for i in range(len(segments)):
    filtered_segment = np.convolve(segments[i], b, mode='full')    # full模式用于保留所有卷积结果 N = L + M -1
    filtered_signals[i*segment_len:i*segment_len+len(filtered_segment)] += filtered_segment # 叠加过程

filtered_signals = filtered_signals#[:target_length] # 保留和原信号等长的部分

# Filter Freq Response
w, h = signal.freqz(b, 1, fs=fs)
plt.figure()
plt.plot(w, abs(h))
plt.title('Filter Freq Response')
plt.grid()
plt.xlabel('f[Hz]')
plt.ylabel('Amplitude')

# Signal Before filtering & Spectrum
n = len(x)
freq = np.fft.fftfreq(n, 1/fs)
y = np.fft.fft(x)

plt.figure()
plt.subplot(221)
plt.plot(t[:500], x[:500])
plt.title('Original Signal')
plt.subplot(222)
plt.plot(freq[:n//2], np.abs(y[:n//2]*2/n)) # 标幺，绘制前一半
plt.title('Spectrum of Orginal Signal')
plt.grid()

# Signal After filtering & Spectrum
n = len(x)
y = np.fft.fft(filtered_signals)

plt.subplot(223)
# 1. Normal output
plt.plot(t[:500], filtered_signals[:500])
plt.title('Filtered Signal')
plt.subplot(224)
plt.plot(freq[:n//2], np.abs(y[:n//2]*2/n)) # 标幺，绘制前一半
plt.title('Spectrum of Filtered Signal')
plt.grid()

plt.tight_layout()

plt.show()

# 2. End of convolution without end cutoff: for test purpose
plt.figure()
temp = t[8000:8850]
s = filtered_signals[8079:8929]
plt.plot(temp, s)
plt.title('Filtered Signal')
plt.grid()
plt.show()