Python:只有 numpy 的低通滤波器
Python: Lowpass Filter with only numpy
我需要在 Python 中实现低通滤波器,但我唯一可以使用的模块是 numpy(不是 scipy)。我尝试在信号上使用 np.fft.fft()
,然后将所有高于截止频率的频率设置为 0,然后使用 np.fft.ifft()
。
但是这不起作用,我根本不知道如何应用过滤器。
编辑:
将 np.abs()
更改为 np.real()
后,结果几乎是正确的。
但是在频谱图中,振幅比原始和过滤后的参考值小(相差 6dB)。所以看起来并不完全正确。有什么想法可以解决这个问题吗?
我的低通函数应该采用以下参数:
signal: audio signal to be filtered
cutoff_freq: cout off frequency in Hz above which to cut off frequencies
sampling_rate: sampling rate in samples/second
应返回过滤后的信号。
我目前的功能
def low_pass_filter(adata: np.ndarray, bandlimit: int = 1000, sampling_rate: int = 44100) -> np.ndarray:
# translate bandlimit from Hz to dataindex according to sampling rate and data size
bandlimit_index = int(bandlimit * adata.size / sampling_rate)
fsig = np.fft.fft(adata)
for i in range(bandlimit_index + 1, len(fsig)):
fsig[i] = 0
adata_filtered = np.fft.ifft(fsig)
return np.real(adata_filtered)
我看到@Cris Luengo 的评论已经将您的解决方案发展到正确的方向。您现在遗漏的最后一件事是,您从 np.fft.fft
获得的频谱由前半部分的正频率分量和后半部分的 'mirrored' 负频率分量组成。
如果您现在将 bandlimit_index
之外的所有分量都设置为零,您将消除这些负频率分量。这解释了 6dB 的信号幅度下降,你消除了一半的信号功率(+ 因为你已经注意到每个真实信号都必须具有共轭对称频谱)。 np.fft.ifft
函数文档 (ifft documentation) 很好地解释了预期的格式。它指出:
“输入的顺序应与 fft 返回的顺序相同,即,”
- a[0]应该包含零频项,
- a[1:n//2] 应该包含 positive-frequency 项,
- a[n//2 + 1:] 应包含 negative-frequency 项,从最负频率开始按递增顺序排列。
这基本上就是您必须保持的对称性。因此,为了保留这些组件,只需将 bandlimit_index + 1 -> (len(fsig) - bandlimit_index)
之间的组件设置为零即可。
def low_pass_filter(adata: np.ndarray, bandlimit: int = 1000, sampling_rate: int = 44100) -> np.ndarray:
# translate bandlimit from Hz to dataindex according to sampling rate and data size
bandlimit_index = int(bandlimit * adata.size / sampling_rate)
fsig = np.fft.fft(adata)
for i in range(bandlimit_index + 1, len(fsig) - bandlimit_index ):
fsig[i] = 0
adata_filtered = np.fft.ifft(fsig)
return np.real(adata_filtered)
或者如果你想要稍微多一些 'pythonic',你也可以像这样将组件设置为零:
fsig[bandlimit_index+1 : -bandlimit_index] = 0
这里有一个 colab 来演示:
https://colab.research.google.com/drive/1RR_9EYlApDMg4jAS2HuJIpSqwg5RLzGW?usp=sharing
实信号的完整 FFT 包含对称的两半。每一半都将包含另一半的反射(纯虚构)复共轭:奈奎斯特之后的一切都不是真实信息。如果您将超过某个频率的所有内容都设置为零,则必须遵守对称性。
这是一个人为设计的信号,以 1kHz 采样,其 FFT:
import numpy as np
from matplotlib import pyplot as plt
t = np.linspace(0, 10, 10000, endpoint=False)
y = np.sin(2 * np.pi * 2 * t) * np.exp(-0.5 * ((t - 5) / 1.5)**2)
f = np.fft.fftfreq(t.size, 0.001)
F = np.fft.fft(y)
plt.figure(constrained_layout=True)
plt.subplot(1, 2, 1)
plt.plot(t, y)
plt.title('Time Domain')
plt.xlabel('time (s)')
plt.subplot(1, 2, 2)
plt.plot(np.fft.fftshift(f), np.fft.fftshift(F.imag), label='imag')
plt.xlim([-3, 3])
plt.title('Frequency Domain')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Imginatry Component')
频率轴是这样的:
>>> f
array([ 0. , 0.1, 0.2, ..., -0.3, -0.2, -0.1])
请注意,除了直流分量 (bin 0) 之外,该轴关于中点是对称的,最高 (奈奎斯特) 频率位于中间。这就是为什么我调用 fftshift
来绘制绘图的原因:它重新排列数组以从最小到最大。
您很可能不需要将输入限制为整数。小数 band_limit
是完全可以接受的。请记住,要将频率转换为索引,您需要乘以大小和采样频率(除以 ram,而不是除以:
def low_pass_filter(data, band_limit, sampling_rate):
cutoff_index = int(band_limit * data.size / sampling_rate)
F = np.fft.fft(data)
F[cutoff_index + 1 : -cutoff_index] = 0
return np.fft.ifft(F).real
您仍然需要 return 实部分量,因为 FFT 在最低的几位中总会有一些虚部舍入误差。
这是截断高于 2Hz 的样本信号图:
y2 = low_pass_filter(y, 2, 1000)
f2 = np.fft.fftfreq(t.size, 0.001)
F2 = np.fft.fft(y2)
plt.figure(constrained_layout=True)
plt.subplot(1, 2, 1)
plt.plot(t, y2)
plt.title('Time Domain')
plt.xlabel('time (s)')
plt.subplot(1, 2, 2)
plt.plot(np.fft.fftshift(f2), np.fft.fftshift(F2.imag), label='imag')
plt.xlim([-3, 3])
plt.title('Frequency Domain')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Imginatry Component')
还记得我们说过纯实数(或纯复数)信号的 FFT 只有一半包含 non-redundant 信息吗? Numpy 尊重这一点,并提供 np.fft.rfft
, np.fft.irfft
, np.fft.rfftfreq
来处理 real-values 信号。您可以使用它来编写一个更简单版本的过滤器,因为不再有对称约束。
def low_pass_filter(data, band_limit, sampling_rate):
cutoff_index = int(band_limit * data.size / sampling_rate)
F = np.fft.rfft(data)
F[cutoff_index + 1:] = 0
return np.fft.irfft(F, n=data.size).real
唯一需要注意的是,我们必须显式地将 n
传递给 irfft
,否则无论输入大小如何,输出大小将始终为偶数。
我需要在 Python 中实现低通滤波器,但我唯一可以使用的模块是 numpy(不是 scipy)。我尝试在信号上使用 np.fft.fft()
,然后将所有高于截止频率的频率设置为 0,然后使用 np.fft.ifft()
。
但是这不起作用,我根本不知道如何应用过滤器。
编辑:
将 np.abs()
更改为 np.real()
后,结果几乎是正确的。
但是在频谱图中,振幅比原始和过滤后的参考值小(相差 6dB)。所以看起来并不完全正确。有什么想法可以解决这个问题吗?
我的低通函数应该采用以下参数:
signal: audio signal to be filtered
cutoff_freq: cout off frequency in Hz above which to cut off frequencies
sampling_rate: sampling rate in samples/second
应返回过滤后的信号。
我目前的功能
def low_pass_filter(adata: np.ndarray, bandlimit: int = 1000, sampling_rate: int = 44100) -> np.ndarray:
# translate bandlimit from Hz to dataindex according to sampling rate and data size
bandlimit_index = int(bandlimit * adata.size / sampling_rate)
fsig = np.fft.fft(adata)
for i in range(bandlimit_index + 1, len(fsig)):
fsig[i] = 0
adata_filtered = np.fft.ifft(fsig)
return np.real(adata_filtered)
我看到@Cris Luengo 的评论已经将您的解决方案发展到正确的方向。您现在遗漏的最后一件事是,您从 np.fft.fft
获得的频谱由前半部分的正频率分量和后半部分的 'mirrored' 负频率分量组成。
如果您现在将 bandlimit_index
之外的所有分量都设置为零,您将消除这些负频率分量。这解释了 6dB 的信号幅度下降,你消除了一半的信号功率(+ 因为你已经注意到每个真实信号都必须具有共轭对称频谱)。 np.fft.ifft
函数文档 (ifft documentation) 很好地解释了预期的格式。它指出:
“输入的顺序应与 fft 返回的顺序相同,即,”
- a[0]应该包含零频项,
- a[1:n//2] 应该包含 positive-frequency 项,
- a[n//2 + 1:] 应包含 negative-frequency 项,从最负频率开始按递增顺序排列。
这基本上就是您必须保持的对称性。因此,为了保留这些组件,只需将 bandlimit_index + 1 -> (len(fsig) - bandlimit_index)
之间的组件设置为零即可。
def low_pass_filter(adata: np.ndarray, bandlimit: int = 1000, sampling_rate: int = 44100) -> np.ndarray:
# translate bandlimit from Hz to dataindex according to sampling rate and data size
bandlimit_index = int(bandlimit * adata.size / sampling_rate)
fsig = np.fft.fft(adata)
for i in range(bandlimit_index + 1, len(fsig) - bandlimit_index ):
fsig[i] = 0
adata_filtered = np.fft.ifft(fsig)
return np.real(adata_filtered)
或者如果你想要稍微多一些 'pythonic',你也可以像这样将组件设置为零:
fsig[bandlimit_index+1 : -bandlimit_index] = 0
这里有一个 colab 来演示: https://colab.research.google.com/drive/1RR_9EYlApDMg4jAS2HuJIpSqwg5RLzGW?usp=sharing
实信号的完整 FFT 包含对称的两半。每一半都将包含另一半的反射(纯虚构)复共轭:奈奎斯特之后的一切都不是真实信息。如果您将超过某个频率的所有内容都设置为零,则必须遵守对称性。
这是一个人为设计的信号,以 1kHz 采样,其 FFT:
import numpy as np
from matplotlib import pyplot as plt
t = np.linspace(0, 10, 10000, endpoint=False)
y = np.sin(2 * np.pi * 2 * t) * np.exp(-0.5 * ((t - 5) / 1.5)**2)
f = np.fft.fftfreq(t.size, 0.001)
F = np.fft.fft(y)
plt.figure(constrained_layout=True)
plt.subplot(1, 2, 1)
plt.plot(t, y)
plt.title('Time Domain')
plt.xlabel('time (s)')
plt.subplot(1, 2, 2)
plt.plot(np.fft.fftshift(f), np.fft.fftshift(F.imag), label='imag')
plt.xlim([-3, 3])
plt.title('Frequency Domain')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Imginatry Component')
频率轴是这样的:
>>> f
array([ 0. , 0.1, 0.2, ..., -0.3, -0.2, -0.1])
请注意,除了直流分量 (bin 0) 之外,该轴关于中点是对称的,最高 (奈奎斯特) 频率位于中间。这就是为什么我调用 fftshift
来绘制绘图的原因:它重新排列数组以从最小到最大。
您很可能不需要将输入限制为整数。小数 band_limit
是完全可以接受的。请记住,要将频率转换为索引,您需要乘以大小和采样频率(除以 ram,而不是除以:
def low_pass_filter(data, band_limit, sampling_rate):
cutoff_index = int(band_limit * data.size / sampling_rate)
F = np.fft.fft(data)
F[cutoff_index + 1 : -cutoff_index] = 0
return np.fft.ifft(F).real
您仍然需要 return 实部分量,因为 FFT 在最低的几位中总会有一些虚部舍入误差。
这是截断高于 2Hz 的样本信号图:
y2 = low_pass_filter(y, 2, 1000)
f2 = np.fft.fftfreq(t.size, 0.001)
F2 = np.fft.fft(y2)
plt.figure(constrained_layout=True)
plt.subplot(1, 2, 1)
plt.plot(t, y2)
plt.title('Time Domain')
plt.xlabel('time (s)')
plt.subplot(1, 2, 2)
plt.plot(np.fft.fftshift(f2), np.fft.fftshift(F2.imag), label='imag')
plt.xlim([-3, 3])
plt.title('Frequency Domain')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Imginatry Component')
还记得我们说过纯实数(或纯复数)信号的 FFT 只有一半包含 non-redundant 信息吗? Numpy 尊重这一点,并提供 np.fft.rfft
, np.fft.irfft
, np.fft.rfftfreq
来处理 real-values 信号。您可以使用它来编写一个更简单版本的过滤器,因为不再有对称约束。
def low_pass_filter(data, band_limit, sampling_rate):
cutoff_index = int(band_limit * data.size / sampling_rate)
F = np.fft.rfft(data)
F[cutoff_index + 1:] = 0
return np.fft.irfft(F, n=data.size).real
唯一需要注意的是,我们必须显式地将 n
传递给 irfft
,否则无论输入大小如何,输出大小将始终为偶数。