为什么在时域和频域中执行卷积结果具有不同的长度?

发布于 2024-09-25 16:50:29 字数 659 浏览 0 评论 0原文

我不是 DSP 专家,但我知道有两种方法可以将离散时域滤波器应用于离散时域波形。第一种方法是在时域中对它们进行卷积,第二种方法是对两者进行 FFT,将两个复频谱相乘,然后对结果进行 IFFT。这些方法的一个关键区别是第二种方法采用循环卷积。

例如,如果滤波器和波形都是 N 点长,则第一种方法(即卷积)会产生 N+N-1 点长的结果,其中此响应的前半部分是滤波器填充,第二部分是滤波器填充。一半是过滤器排空。为了获得稳态响应,滤波器的点数需要少于要滤波的波形的点数。

继续使用第二种方法的示例,并假设离散时域波形数据都是实数(不复杂),滤波器的 FFT 和波形都产生 N 点长的 FFT。将两个频谱 IFFT 的结果相乘会产生一个同样为 N 点长的时域结果。这里,滤波器充满和排空的响应在时域中相互重叠,并且没有稳态响应。这就是循环卷积的效果。为了避免这种情况,滤波器尺寸通常会小于波形尺寸,并且两者都将被补零,以允许频率卷积空间在两个频谱乘积的 IFFT 之后及时扩展。

我的问题是,我经常在文献中看到知名专家/公司的工作,他们有一个离散(真实)时域波形(N 点),他们对其进行 FFT,乘以一些滤波器(也是 N 点),对结果进行IFFT以供后续处理。我天真的想法是这个结果应该不包含稳态响应,因此应该包含过滤器填充/清空的伪影,这会导致解释结果数据时出现错误,但我一定遗漏了一些东西。在什么情况下这是一个有效的方法?

任何见解将不胜感激

I'm not a DSP expert, but I understand that there are two ways that I can apply a discrete time-domain filter to a discrete time-domain waveform. The first is to convolve them in the time domain, and the second is to take the FFT of both, multiply both complex spectrums, and take IFFT of the result. One key difference in these methods is the second approach is subject to circular convolution.

As an example, if the filter and waveforms are both N points long, the first approach (i.e. convolution) produces a result that is N+N-1 points long, where the first half of this response is the filter filling up and the 2nd half is the filter emptying. To get a steady-state response, the filter needs to have fewer points than the waveform to be filtered.

Continuing this example with the second approach, and assuming the discrete time-domain waveform data is all real (not complex), the FFT of the filter and the waveform both produce FFTs of N points long. Multiplying both spectrums IFFT'ing the result produces a time-domain result also N points long. Here the response where the filter fills up and empties overlap each other in the time domain, and there's no steady state response. This is the effect of circular convolution. To avoid this, typically the filter size would be smaller than the waveform size and both would be zero-padded to allow space for the frequency convolution to expand in time after IFFT of the product of the two spectrums.

My question is, I often see work in the literature from well-established experts/companies where they have a discrete (real) time-domain waveform (N points), they FFT it, multiply it by some filter (also N points), and IFFT the result for subsequent processing. My naive thinking is this result should contain no steady-state response and thus should contain artifacts from the filter filling/emptying that would lead to errors in interpreting the resulting data, but I must be missing something. Under what circumstances can this be a valid approach?

Any insight would be greatly appreciated

如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。

评论(4

长亭外,古道边 2024-10-02 16:50:30

尽管假设矩形数据窗口在 FFT 孔径宽度处是周期性的(这是对循环卷积在没有足够零填充的情况下所做的事情的一种解释)会产生伪影,但差异可能大到足以淹没数据分析。问题。

Although there will be artifacts from assuming that a rectangular window of data is periodic at the FFT aperture width, which is one interpretation of what circular convolution does without sufficient zero padding, the differences may or may not be large enough to swamp the data analysis in question.

肤浅与狂妄 2024-10-02 16:50:29

基本问题不在于零填充与假设的周期性,而是傅里叶分析将信号分解为正弦波,在最基本的层面上,假设其范围是无限的。 两种方法都是正确的,因为使用完整 FFT 的 IFFT 将返回精确的输入波形,而两种方法都是错误的,因为使用少于全频谱的频谱可能会导致影响在边缘(通常延伸几个波长)。唯一的区别在于你假设填充无穷大的其余部分的细节,而不在于你是否做出假设。

回到你的第一段:通常,在 DSP 中,我遇到的 FFT 最大问题是它们是非因果的,因此我通常更喜欢留在时域,使用例如 FIR 和 IIR 滤波器。

更新:

在问题陈述中,OP正确指出了使用 FFT 过滤信号时可能出现的一些问题,例如边缘效应,在进行以下卷积时可能会特别成问题:长度(在时域)与采样波形相当。值得注意的是,并非所有过滤都是使用 FFT 完成的,并且在 OP 引用的论文中,他们没有使用 FFT 滤波器,并且使用他们的方法不会出现 FFT 滤波器实现中出现的问题。

例如,考虑一个使用两种不同实现方式对 128 个采样点进行简单平均的滤波器。

FFT:在 FFT/卷积方法中,我们将拥有 256 个点的样本,并将其与前半部分恒定并在后半部分变为零的 wfm 进行卷积。这里的问题是(即使这个系统已经运行了几个周期),什么决定了结果的第一个点的值? FFT 假设 wfm 是圆形的(即无限周期性),因此:结果的第一个点由 wfm 的最后 127 个(即未来)样本确定(跳过wfm),或者如果进行零填充则增加 127 个零。两者都不正确。

FIR:另一种方法是使用 FIR 滤波器实现平均值。例如,这里可以使用 128 个寄存器 FIFO 队列中的值的平均值。也就是说,当每个样本点进入时,1)将其放入队列中,2)将最旧的项目出列,3)对队列中剩余的所有 128 个项目进行平均;这是该样本点的结果。这种方法连续运行,一次处理一个点,并在每个样本后返回过滤后的结果,并且不会出现 FFT 应用于有限样本块时出现的任何问题。每个结果只是当前样本和之前 127 个样本的平均值。

OP 引用的论文采用的方法与 FIR 滤波器比 FFT 滤波器更相似(请注意,论文中的滤波器更复杂,整篇论文基本上都是对该滤波器的分析。)例如,参见,这本免费书籍,它描述了如何分析和应用不同的滤波器,还请注意拉普拉斯方法来分析FIR 和 IIR 滤波器与引用论文中的内容非常相似。

The basic problem is not about zero padding vs the assumed periodicity, but that Fourier analysis decomposes the signal into sine waves which, at the most basic level, are assumed to be infinite in extent. Both approaches are correct in that the IFFT using the full FFT will return the exact input waveform, and both approaches are incorrect in that using less than the full spectrum can lead to effects at the edges (that usually extend a few wavelengths). The only difference is in the details of what you assume fills in the rest of infinity, not in whether you are making an assumption.

Back to your first paragraph: Usually, in DSP, the biggest problem I run into with FFTs is that they are non-causal, and for this reason I often prefer to stay in the time domain, using, for example, FIR and IIR filters.

Update:

In the question statement, the OP correctly points out some of the problems that can arise when using FFTs to filter signals, for example, edge effects, that can be particularly problematic when doing a convolution that is comparable in the length (in the time domain) to the sampled waveform. It's important to note though that not all filtering is done using FFTs, and in the paper cited by the OP, they are not using FFT filters, and the problems that would arise with an FFT filter implementation do not arise using their approach.

Consider, for example, a filter that implements a simple average over 128 sample points, using two different implementations.

FFT: In the FFT/convolution approach one would have a sample of, say, 256, points and convolve this with a wfm that is constant for the first half and goes to zero in the second half. The question here is (even after this system has run a few cycles), what determines the value of the first point of the result? The FFT assumes that the wfm is circular (i.e. infinitely periodic) so either: the first point of the result is determined by the last 127 (i.e. future) samples of the wfm (skipping over the middle of the wfm), or by 127 zeros if you zero-pad. Neither is correct.

FIR: Another approach is to implement the average with an FIR filter. For example, here one could use the average of the values in a 128 register FIFO queue. That is, as each sample point comes in, 1) put it in the queue, 2) dequeue the oldest item, 3) average all of the 128 items remaining in the queue; and this is your result for this sample point. This approach runs continuously, handling one point at a time, and returning the filtered result after each sample, and has none of the problems that occur from the FFT as it's applied to finite sample chunks. Each result is just the average of the current sample and the 127 samples that came before it.

The paper that OP cites takes an approach much more similar to the FIR filter than to the FFT filter (note though that the filter in the paper is more complicated, and the whole paper is basically an analysis of this filter.) See, for example, this free book which describes how to analyze and apply different filters, and note also that the Laplace approach to analysis of the FIR and IIR filters is quite similar what what's found in the cited paper.

面如桃花 2024-10-02 16:50:29

下面是 DFT(循环卷积)与线性卷积的无零填充卷积示例。这是长度 M=32 的序列与长度 L=128 的序列的卷积(使用 Numpy/Matplotlib):

f = rand(32); g = rand(128)
h1 = convolve(f, g)
h2 = real(ifft(fft(f, 128)*fft(g)))
plot(h1); plot(h2,'r')
grid()

alt text
前 M-1 点不同,由于没有用零填充,因此少了 M-1 点。如果您正在进行块卷积,这些差异是一个问题,但诸如重叠和保存或重叠和相加等技术可以用来克服这个问题。否则,如果您只是计算一次性过滤操作,则有效结果将从索引 M-1 开始,到索引 L-1 结束,长度为 L-M+1。

至于引用的论文,我在附录 A 中查看了他们的 MATLAB 代码。我认为他们在将 Hfinal 传递函数应用于负频率而不首先将其共轭时犯了一个错误。否则,您可以在他们的图表中看到时钟抖动是周期性信号,因此使用循环卷积对于稳态分析来说是很好的。

编辑:关于共轭传递函数,PLL 具有实值脉冲响应,并且每个实值信号具有共轭对称频谱。在代码中您可以看到他们只是使用 Hfinal[Ni] 来获取负频率而不采用共轭。我绘制了从 -50 MHz 到 50 MHz 的传递函数:

N = 150000                    # number of samples. Need >50k to get a good spectrum. 
res = 100e6/N                 # resolution of single freq point  
f = res * arange(-N/2, N/2)   # set the frequency sweep [-50MHz,50MHz), N points
s = 2j*pi*f                   # set the xfer function to complex radians 

f1 = 22e6       # define 3dB corner frequency for H1 
zeta1 = 0.54    # define peaking for H1 
f2 = 7e6        # define 3dB corner frequency for H2 
zeta2 = 0.54    # define peaking for H2    
f3 = 1.0e6      # define 3dB corner frequency for H3 

# w1 = natural frequency   
w1 = 2*pi*f1/((1 + 2*zeta1**2 + ((1 + 2*zeta1**2)**2 + 1)**0.5)**0.5)  
# H1 transfer function 
H1 = ((2*zeta1*w1*s + w1**2)/(s**2 + 2*zeta1*w1*s + w1**2))            

# w2 = natural frequency 
w2 = 2*pi*f2/((1 + 2*zeta2**2 + ((1 + 2*zeta2**2)**2 + 1)**0.5)**0.5)  
# H2 transfer function  
H2 = ((2*zeta2*w2*s + w2**2)/(s**2 + 2*zeta2*w2*s + w2**2))            

w3 = 2*pi*f3        # w3 = 3dB point for a single pole high pass function. 
H3 = s/(s+w3)       # the H3 xfer function is a high pass

Ht = 2*(H1-H2)*H3   # Final transfer based on the difference functions

subplot(311); plot(f, abs(Ht)); ylabel("abs")
subplot(312); plot(f, real(Ht)); ylabel("real")
subplot(313); plot(f, imag(Ht)); ylabel("imag")

alt text

如您所见,实际组件具有偶对称性,虚部具有奇对称性。在他们的代码中,他们只计算对数图的正频率(足够合理)。然而,为了计算逆变换,他们通过索引 Hfinal[Ni] 使用正频率的值作为负频率,但忘记将其共轭。

Here's an example of convolution without zero padding for the DFT (circular convolution) vs linear convolution. This is the convolution of a length M=32 sequence with a length L=128 sequence (using Numpy/Matplotlib):

f = rand(32); g = rand(128)
h1 = convolve(f, g)
h2 = real(ifft(fft(f, 128)*fft(g)))
plot(h1); plot(h2,'r')
grid()

alt text
The first M-1 points are different, and it's short by M-1 points since it wasn't zero padded. These differences are a problem if you're doing block convolution, but techniques such as overlap and save or overlap and add are used to overcome this problem. Otherwise if you're just computing a one-off filtering operation, the valid result will start at index M-1 and end at index L-1, with a length of L-M+1.

As to the paper cited, I looked at their MATLAB code in appendix A. I think they made a mistake in applying the Hfinal transfer function to the negative frequencies without first conjugating it. Otherwise, you can see in their graphs that the clock jitter is a periodic signal, so using circular convolution is fine for a steady-state analysis.

Edit: Regarding conjugating the transfer function, the PLLs have a real-valued impulse response, and every real-valued signal has a conjugate symmetric spectrum. In the code you can see that they're just using Hfinal[N-i] to get the negative frequencies without taking the conjugate. I've plotted their transfer function from -50 MHz to 50 MHz:

N = 150000                    # number of samples. Need >50k to get a good spectrum. 
res = 100e6/N                 # resolution of single freq point  
f = res * arange(-N/2, N/2)   # set the frequency sweep [-50MHz,50MHz), N points
s = 2j*pi*f                   # set the xfer function to complex radians 

f1 = 22e6       # define 3dB corner frequency for H1 
zeta1 = 0.54    # define peaking for H1 
f2 = 7e6        # define 3dB corner frequency for H2 
zeta2 = 0.54    # define peaking for H2    
f3 = 1.0e6      # define 3dB corner frequency for H3 

# w1 = natural frequency   
w1 = 2*pi*f1/((1 + 2*zeta1**2 + ((1 + 2*zeta1**2)**2 + 1)**0.5)**0.5)  
# H1 transfer function 
H1 = ((2*zeta1*w1*s + w1**2)/(s**2 + 2*zeta1*w1*s + w1**2))            

# w2 = natural frequency 
w2 = 2*pi*f2/((1 + 2*zeta2**2 + ((1 + 2*zeta2**2)**2 + 1)**0.5)**0.5)  
# H2 transfer function  
H2 = ((2*zeta2*w2*s + w2**2)/(s**2 + 2*zeta2*w2*s + w2**2))            

w3 = 2*pi*f3        # w3 = 3dB point for a single pole high pass function. 
H3 = s/(s+w3)       # the H3 xfer function is a high pass

Ht = 2*(H1-H2)*H3   # Final transfer based on the difference functions

subplot(311); plot(f, abs(Ht)); ylabel("abs")
subplot(312); plot(f, real(Ht)); ylabel("real")
subplot(313); plot(f, imag(Ht)); ylabel("imag")

alt text

As you can see, the real component has even symmetry and the imaginary component has odd symmetry. In their code they only calculated the positive frequencies for a loglog plot (reasonable enough). However, for calculating the inverse transform they used the values for the positive frequencies for the negative frequencies by indexing Hfinal[N-i] but forgot to conjugate it.

鯉魚旗 2024-10-02 16:50:29

我可以阐明为什么在应用 FFT 之前应用“加窗”的原因。

正如已经指出的,FFT 假设我们有一个无限信号。当我们在有限时间 T 内采样时,这在数学上相当于将信号乘以矩形函数。

时域中的乘法变为频域中的卷积。矩形的频率响应是同步函数,即sin(x)/x。分子中的 x 是最重要的,因为它会减少 O(1/N)。

如果您的频率分量正好是 1/T 的倍数,这并不重要,因为同步函数在除频率为 1 之外的所有点上均为零。

但是,如果您有一个落在 2 点之间的正弦,您将看到同步在频率点上采样的函数。它看起来像是同步函数的放大版本,并且由卷积引起的“幽灵”信号以 1/N 或 6dB/倍频程衰减。如果您的信号高于本底噪声 60db,您将看不到主信号左右 1000 个频率的噪声,它将被同步功能的“裙子”淹没。

如果您使用不同的时间窗口,您会得到不同的频率响应,例如余弦随 1/x^2 减弱,有针对不同测量的专用窗口。汉宁窗通常用作通用窗。

关键是,在不应用任何“窗口函数”时使用的矩形窗口会比精心选择的窗口产生更糟糕的伪影。即,通过“扭曲”时间样本,我们在频域中获得更好的图像,更接近“现实”,或者更确切地说,我们期望并希望看到的“现实”。

I can shed some light to the reason why "windowing" is applied before FFT is applied.

As already pointed out the FFT assumes that we have a infinite signal. When we take a sample over a finite time T this is mathematically the equivalent of multiplying the signal with a rectangular function.

Multiplying in the time domain becomes convolution in the frequency domain. The frequency response of a rectangle is the sync function i.e. sin(x)/x. The x in the numerator is the kicker, because it dies down O(1/N).

If you have frequency components which are exactly multiples of 1/T this does not matter as the sync function is zero in all points except that frequency where it is 1.

However if you have a sine which fall between 2 points you will see the sync function sampled on the frequency point. It lloks like a magnified version of the sync function and the 'ghost' signals caused by the convolution die down with 1/N or 6dB/octave. If you have a signal 60db above the noise floor, you will not see the noise for 1000 frequencies left and right from your main signal, it will be swamped by the "skirts" of the sync function.

If you use a different time window you get a different frequency response, a cosine for example dies down with 1/x^2, there are specialized windows for different measurements. The Hanning window is often used as a general purpose window.

The point is that the rectangular window used when not applying any "windowing function" creates far worse artefacts than a well chosen window. i.e by "distorting" the time samples we get a much better picture in the frequency domain which closer resembles "reality", or rather the "reality" we expect and want to see.

~没有更多了~
我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
原文