滤波 - 频域乘法

信息处理 fft 过滤器 频谱 Python
2021-12-25 01:02:42

我正在尝试创建一个简单的低通滤波器,但是在查看一个简单的巴特沃斯滤波器的频率响应时,我得到了一个令人惊讶的结果。

我从另一篇文章中复制了下面的大部分示例。我在脚本底部添加了一些代码,以将输入和输出频谱与滤波器的频率响应进行比较。我希望输出频谱 B 应该是输入频谱 A 和频率响应 H 的乘积:

B=HA

然而,下图显示滤波器实际上增加了一些低频分量 - 看看红线是如何在 4 Hz 附近的绿色上方。

谁能解释这是为什么?

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt
from scipy.fftpack import fft as fft
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y


# Filter requirements.
order = 6
fs = 30.0       # sample rate, Hz
cutoff = 3.667  # desired cutoff frequency of the filter, Hz

# Get the filter coefficients so we can check its frequency response.
b, a = butter_lowpass(cutoff, fs, order)

# Plot the frequency response.
w, h = freqz(b, a, worN=8000)
plt.subplot(2, 1, 1)
plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b')
plt.plot(cutoff, 0.5*np.sqrt(2), 'ko')
plt.axvline(cutoff, color='k')
plt.xlim(0, 0.5*fs)
plt.title("Lowpass Filter Frequency Response")
plt.xlabel('Frequency [Hz]')
plt.grid()


# Demonstrate the use of the filter.
# First make some data to be filtered.
T = 5.0         # seconds
n = int(T * fs) # total number of samples
t = np.linspace(0, T, n, endpoint=False)
# "Noisy" data.  We want to recover the 1.2 Hz signal from this.
data = np.sin(1.2*2*np.pi*t) + 1.5*np.cos(9*2*np.pi*t) + 0.5*np.sin(12.0*2*np.pi*t)

# Filter the data, and plot both the original and filtered signals.
y = butter_lowpass_filter(data, cutoff, fs, order)

plt.subplot(2, 1, 2)
plt.plot(t, data, 'b-', label='data')
plt.plot(t, y, 'g-', linewidth=2, label='filtered data')
plt.xlabel('Time [sec]')
plt.grid()
plt.legend()

plt.subplots_adjust(hspace=0.35)
plt.show()

def calculateFFT(time,signal):
    N=len(signal)      
    df=1/((time[-1]-time[0]))
    frequencies=[i*df for i in range(int(N/2.0))]
    fftValues = [2.0/N*abs(i) for i in fft(signal,N)[0:N/2.0] ]
    return frequencies,fftValues

plt.subplot(2, 1, 1)
originalfreqs,originalFFT=calculateFFT(t,data)
plt.plot(originalfreqs,originalFFT,"g",label="original")

filteredfreqs,filteredFFT=calculateFFT(t,y)
plt.plot(filteredfreqs,filteredFFT,"r",label="filtered")
plt.legend()

在此处输入图像描述

2个回答

要扩展@Maximilian Matthé答案,您可以通过对 中使用的输入进行零填充来可视化频谱泄漏效应(频域中 sinc 函数的卷积)calculateFFT例如,以下函数将输入零填充到比原始输入长 $k$ 倍的长度(在本例中 $k=4$):k times longer than the original input (where in this instance k=4):

def calculateFFT(time,signal):
    k=4
    N=k*len(signal)      
    df=1/(k*(time[-1]-time[0]))
    frequencies=[i*df for i in range(int(N/2.0))]
    fftValues = [k*2.0/N*abs(i) for i in fft(signal,N)[0:int(N/2.0)] ]
    return frequencies,fftValues

绘制结果,您应该看到低频分量实际上重叠(这表明滤波实际上并没有增加那里的信号):

在此处输入图像描述

那么,尖刺周围的那些涟漪是从哪里来的呢?事实证明它们一直都在那里,但是在计算 FFT 时,您会在一组离散的频率值处获得频谱值​​,而这些波纹恰好在这些频率处通过零。如果您选择了一个稍微不同的信号频率,而不是 FFT 频率分辨率的精确倍数(在您的情况下为 0.2Hz),例如 1.25Hz 而不是 1.2Hz,则频谱的采样看起来会完全不同(由于到波纹振荡中不同点的频率采样):

在此处输入图像描述

我很高兴地支持这个问题,因为这是一个好问题应该是什么样子,也是我希望学生如何验证他们对所学内容的理解。热衷于了解您以后想要使用的东西的背景总是好的。

您遇到的问题如下:原则上您是对的,卷积定理意味着一个域中的卷积是另一个域中的乘法。因此你有

F{x(t)h(t)}=X(f)H(f)
where F denotes the Fourier Transform.

那么,您的系统中的 $x(t)$ 和 $h(t)$ 是什么?$x(t)$ 是您的输入信号(蓝色),它是三个正弦的总和,乘以一个矩形窗口。它隐含地乘以一个矩形,因为您的时间(在模拟中,它不能)从 $-\infty$ 到 $\infty$。因此,$X(f)$ 的频谱是 sinc 函数(来自矩形窗口)与三个狄拉克(在正弦频率处)的卷积。显然,这不是一个纯粹的离散谱。x(t) and h(t) in your system? x(t) is your input signal (in blue), which is the sum of three sines, multiplied by a rectangular window. It is implicitely multiplied by a rectangular, because your time does not (in simulation, it cannot) reach from to . So, the spectrum of X(f) is the convolution of a sinc function (from the rectangular window) with three Diracs (at the frequency of the sines). Clearly, this is not a purely discrete spectrum.

什么是 $h(t)$?它是巴特沃斯滤波器的脉冲响应。Butterworth 滤波器具有无限长的脉冲响应,因此您的卷积乘积 $x(t)*h(t)$ 也是无限宽的。因此,原则上您不能应用有限(即离散)傅里叶变换(fft)并期望它类似于连续情况。h(t)? It is the impulse response of the butterworth filter. The Butterworth filter has an infinitely long impulse response, hence your convolution product x(t)h(t) is also infinitely wide. So, in principle you cannot apply a finite (i.e. discrete) Fourier Transform (fft) and expect it be a similar to the continuous case.

所以,你看到的是合理的。您可以尝试对输入信号进行零填充,以便获得信号中的脉冲响应(主要部分)。但是,即使那样,您也会看到频率,因为 $X(f)$ 实际上不是离散的。X(f) is not discrete actually.