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

ggk*_*ath 10 signal-processing fft fftw

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

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

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

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

任何见解将不胜感激

tom*_*m10 8

基本问题不是关于零填充与假定的周期性,而是傅里叶分析将信号分解为正弦波,在最基本的水平上,假设在范围内是无限的. 这两种方法都是正确的使用完整的FFT将返回精确的输入波形的IFFT,并且这两种方法是不正确的使用比全谱少可能会导致在边缘效应在(通常延长几个波长).唯一的区别在于你假设填充剩余无穷大的细节,而不是你是否做出假设.

回到你的第一段:通常情况下,在DSP,最大的问题,我跑与FFT的是,他们都是非因果关系,基于这个原因,我经常喜欢留在时域,使用,例如,FIR和IIR滤波器.

更新:

在问题语句中,OP正确地指出某些使用的FFT时进行过滤的信号,可能出现例如存在的问题,边缘效应,做卷积是在长度相当(在时域中的时候,可以是特别有问题的)到采样波形.要注意,虽然并非所有的过滤是使用FFT的完成是非常重要的,并且由OP引用的论文中,他们没有使用FFT过滤器,并使用他们的方法不会出现将与FFT滤波器的实现出现的问题.

例如,考虑使用两种不同的实现方式实现简单平均超过128个采样点的过滤器.

FFT:在FFT /卷积办法一个人必须的,比方说,256分和卷积这个样品与WFM是上半年不变,并在下半年趋于零.这里的问题是(即使在这个系统运行了几个周期之后),是什么决定了结果的第一个点的值?在FFT假定WFM是圆形的(即无限周期的),以便任一:结果的第一点由所确定的最后 127(即,未来)的WFM的样品(跳过的WFM的中间),或者通过127个零如果你零填充.两者都不正确.

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

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


Ery*_*Sun 7

下面是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()
Run Code Online (Sandbox Code Playgroud)

替代文字 第一个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")
Run Code Online (Sandbox Code Playgroud)

替代文字

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