2024-06-12 08:03:50 发布
网友
我想用python来得到一个带通滤波器,它有一个128点的Hamming窗口,截止频率为0.7-4Hz。我从图像中获取信号样本。(1个样本=1个图像)。fps经常变化。
在python中如何实现这一点?我读到这个:http://mpastell.com/2010/01/18/fir-with-scipy/但是我发现firwin相当混乱。这个变量fps怎么能做到?
可以使用函数^{}或^{}创建带通FIR滤波器。也可以使用^{}设计FIR滤波器
下面的代码为创建带通FIR滤波器提供了一些方便的包装。它使用这些来创建与问题中请求的数字对应的带通滤波器。这假设采样是均匀的。如果采样不均匀,则不适合使用FIR滤波器。
from scipy.signal import firwin, remez, kaiser_atten, kaiser_beta # Several flavors of bandpass FIR filters. def bandpass_firwin(ntaps, lowcut, highcut, fs, window='hamming'): nyq = 0.5 * fs taps = firwin(ntaps, [lowcut, highcut], nyq=nyq, pass_zero=False, window=window, scale=False) return taps def bandpass_kaiser(ntaps, lowcut, highcut, fs, width): nyq = 0.5 * fs atten = kaiser_atten(ntaps, width / nyq) beta = kaiser_beta(atten) taps = firwin(ntaps, [lowcut, highcut], nyq=nyq, pass_zero=False, window=('kaiser', beta), scale=False) return taps def bandpass_remez(ntaps, lowcut, highcut, fs, width): delta = 0.5 * width edges = [0, lowcut - delta, lowcut + delta, highcut - delta, highcut + delta, 0.5*fs] taps = remez(ntaps, edges, [0, 1, 0], Hz=fs) return taps if __name__ == "__main__": import numpy as np import matplotlib.pyplot as plt from scipy.signal import freqz # Sample rate and desired cutoff frequencies (in Hz). fs = 63.0 lowcut = 0.7 highcut = 4.0 ntaps = 128 taps_hamming = bandpass_firwin(ntaps, lowcut, highcut, fs=fs) taps_kaiser16 = bandpass_kaiser(ntaps, lowcut, highcut, fs=fs, width=1.6) taps_kaiser10 = bandpass_kaiser(ntaps, lowcut, highcut, fs=fs, width=1.0) remez_width = 1.0 taps_remez = bandpass_remez(ntaps, lowcut, highcut, fs=fs, width=remez_width) # Plot the frequency responses of the filters. plt.figure(1, figsize=(12, 9)) plt.clf() # First plot the desired ideal response as a green(ish) rectangle. rect = plt.Rectangle((lowcut, 0), highcut - lowcut, 1.0, facecolor="#60ff60", alpha=0.2) plt.gca().add_patch(rect) # Plot the frequency response of each filter. w, h = freqz(taps_hamming, 1, worN=2000) plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="Hamming window") w, h = freqz(taps_kaiser16, 1, worN=2000) plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="Kaiser window, width=1.6") w, h = freqz(taps_kaiser10, 1, worN=2000) plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="Kaiser window, width=1.0") w, h = freqz(taps_remez, 1, worN=2000) plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="Remez algorithm, width=%.1f" % remez_width) plt.xlim(0, 8.0) plt.ylim(0, 1.1) plt.grid(True) plt.legend() plt.xlabel('Frequency (Hz)') plt.ylabel('Gain') plt.title('Frequency response of several FIR filters, %d taps' % ntaps) plt.show()
这是脚本生成的情节。当然,在本地运行脚本更有用,因此可以放大细节。
另一个选择是(异步)采样率转换,在过滤之前将数据转换为恒定采样率(如“网格化”)。当然,这只有在你知道采样率的情况下才有效,而且只有在你真的需要过滤数据(而不仅仅是估计光谱)时才有用。
为此,例如,可以逐帧应用scipy.interpolate中的InterpolatedUnivariateSpline,它相当快。
尝试以不一致的采样率过滤数据是非常困难的(不可能?)。所以你要做的是:
创建一个固定采样率的新信号。固定采样率应为最大采样率或更高采样率。要做到这一点,需要设置一个新的“网格”来表示新样本应该去的地方,并从现有数据中插值它们的值。根据需要的精度,存在多种插值方法。线性插值可能不是一个坏的起点,但它取决于你在做什么。如果您不确定,请在https://dsp.stackexchange.com/上询问。
一旦你做到了这一点,你就可以对你的信号应用标准的信号处理方法,因为样本是均匀放置的,比如你链接的帖子中描述的那些。
如有必要,可能需要再次插值以恢复原始采样位置。
如果您只想分析数据,那么您可能对Lomb Periodigram感兴趣。你可以使用Lomb Periodigram,而不是通过频带传递数据然后进行分析,然后只查看相关的频率,或者根据你的需要对结果进行加权。(另请参见数字配方系列。第13.8章被称为“不均匀间隔数据的谱分析”,这似乎是比维基百科页面更友好的介绍)
可以使用函数^{} 或^{} 创建带通FIR滤波器。也可以使用^{} 设计FIR滤波器
下面的代码为创建带通FIR滤波器提供了一些方便的包装。它使用这些来创建与问题中请求的数字对应的带通滤波器。这假设采样是均匀的。如果采样不均匀,则不适合使用FIR滤波器。
这是脚本生成的情节。当然,在本地运行脚本更有用,因此可以放大细节。
另一个选择是(异步)采样率转换,在过滤之前将数据转换为恒定采样率(如“网格化”)。当然,这只有在你知道采样率的情况下才有效,而且只有在你真的需要过滤数据(而不仅仅是估计光谱)时才有用。
为此,例如,可以逐帧应用scipy.interpolate中的InterpolatedUnivariateSpline,它相当快。
尝试以不一致的采样率过滤数据是非常困难的(不可能?)。所以你要做的是:
创建一个固定采样率的新信号。固定采样率应为最大采样率或更高采样率。要做到这一点,需要设置一个新的“网格”来表示新样本应该去的地方,并从现有数据中插值它们的值。根据需要的精度,存在多种插值方法。线性插值可能不是一个坏的起点,但它取决于你在做什么。如果您不确定,请在https://dsp.stackexchange.com/上询问。
一旦你做到了这一点,你就可以对你的信号应用标准的信号处理方法,因为样本是均匀放置的,比如你链接的帖子中描述的那些。
如有必要,可能需要再次插值以恢复原始采样位置。
如果您只想分析数据,那么您可能对Lomb Periodigram感兴趣。你可以使用Lomb Periodigram,而不是通过频带传递数据然后进行分析,然后只查看相关的频率,或者根据你的需要对结果进行加权。(另请参见数字配方系列。第13.8章被称为“不均匀间隔数据的谱分析”,这似乎是比维基百科页面更友好的介绍)
相关问题 更多 >
编程相关推荐