为什么我通过频域中的零填充进行的时域插值是错误的?

信息处理 自由度 插值 零填充
2022-01-30 14:30:37

由于该过程可以应用于任一域以增加另一个域的采样率,因此我尝试在频率空间中应用零填充以在时间空间中恢复“更清晰”的插值信号。为此,我在频谱中较高频率的位置插入零值频率,这是一种常见的做法。

但是,在零填充(红色)之后,我似乎不能很好地恢复原始信号(下面的黑色)。

import numpy as np
import matplotlib.pyplot as plt

# odd dimension for simplicity
n    = 19
npad = 99

x    = np.linspace(0.,4.*np.pi,n)
xpad = np.linspace(0.,4.*np.pi,npad)

f = np.cos(x) + 1j*np.sin(x)

f_fwd = np.fft.fft(f)

f_fwd_pad = np.zeros(npad,dtype=complex)

h = (n-1)//2
f_fwd_pad[0:h+1]   = f_fwd[0:h+1]
f_fwd_pad[npad-h:] = f_fwd[h+1:]

f_interpolated = np.fft.ifft(f_fwd_pad)*npad/n

fig, ax = plt.subplots(1,2)

ax[0].plot(x,np.real(f),linestyle=None,marker='x',color='k')
ax[0].plot(xpad,np.real(f_interpolated),color='r')

ax[1].plot(x,np.imag(f),linestyle=None,marker='x',color='k')
ax[1].plot(xpad,np.imag(f_interpolated),color='r')

在此处输入图像描述

这个结果是预期的吗?是否有一些我缺少的基本理解?

3个回答

上面的答案是正确的。只是为了进一步澄清, usingx = np.linspace(0,10,5)将产生从 0到 10 的 5 个数字(含)

np.linspace(0,10,5)
array([ 0. ,  2.5,  5. ,  7.5, 10. ])

您不想要最后一个数字,因为在您的示例中,最后一个数字是下一个时期的第一个数字。正确的实现是:

periods = 4.*np.pi
x    = np.arange(0., periods, periods/n)
xpad = np.arange(0., periods, periods/npad)

结果与适当的时期

你的问题在这里:

x    = np.linspace(0.,4.*np.pi,n)

周期应该是从 0 到 2pi 的 N+1 个样本,因为 n=0..N。然后取 x(k) for k=0..N-1

目前您的 FFT 不是纯单音,因为正弦曲线在 FFT 周期内没有完美周期。因此用零填充将不是正确的填充。上面的修复将是 FFT 周期内的完美周期,使零填充成为正确的填充。

这是应用了 MG 修复程序的代码。其他一些调整和蓬松的云。这些照片在我看来很好。命令行 Python 2.7

将 numpy 导入为 np
将 matplotlib.pyplot 导入为 plt

# 为简单起见奇数维度
n = 19
npad = 299

持续时间 = 0.85*2.0*np.pi
x = np.arange(0., 持续时间, 持续时间/n)
xpad = np.arange(0., 持续时间, 持续时间/npad)

f = np.cos(x) + 1j*np.sin(x)

f_fwd = np.fft.fft(f)

f_fwd_pad = np.zeros(npad,dtype=complex)

h = (n-1)/2
h = n - 3
f_fwd_pad[0:h+1] = f_fwd[0:h+1]
f_fwd_pad[npad-h:] = f_fwd[nh:]

f_interpolated = np.fft.ifft(f_fwd_pad)*npad/n

无花果,斧头 = plt.subplots(1,2)

ax[0].plot(x,np.real(f),marker='x',color='k')
ax[0].plot(xpad,np.real(f_interpolated),color='r')

ax[1].plot(x,np.imag(f),marker='x',color='k')
ax[1].plot(xpad,np.imag(f_interpolated),color='r')

plt.show()

plt.plot(np.real(f_interpolated),np.imag(f_interpolated))
plt.show()

在此处输入图像描述