在matlab中通过非整数移位来移动向量的元素

Question

我想通过非整数移位来移动向量，线性插值似乎不是很准确，所以我尝试通过以下使用傅立叶变换的代码来使用sinc插值。

function y = fshift(x,s)
% FSHIFT Fractional circular shift
%   Syntax:
%
%       >> y = fshift(x,s)
%
%   FSHIFT circularly shifts the elements of vector x by a (possibly
%   non-integer) number of elements s. FSHIFT works by applying a linear
%   phase in the spectrum domain and is equivalent to CIRCSHIFT for integer
%   values of argument s (to machine precision).


needtr = 0; if size(x,1) == 1; x = x(:); needtr = 1; end;
N = size(x,1); 
r = floor(N/2)+1; f = ((1:N)-r)/(N/2); 
p = exp(-j*s*pi*f)'; 
y = ifft(fft(x).*ifftshift(p)); if isreal(x); y = real(y); end;
if needtr; y = y.'; end;

当我将方波移位整数移位时，代码中没有问题，但当移位为非整数时，输出会出现显着波动 即，

s=[zeros(1,20) ones(1,20) zeros(1,20)];
b=fshift(s,3.5);
stem(b)

如何克服这个问题，还有其他准确的方法吗？

Answer 1

试试这个：

假设您要移动 3.5。找出过采样值是多少（即哪个值会将移位更改为整数 - 在本例中为 2）。

ov = 2;
a = 3.5*ov;
y = downsample(circshift(interp(s,2).', -a),ov);

这在边缘仍然有一些振铃，但比你的 sinc 方法少得多。我不确定这是否是由于吉布斯现象造成的，因为正如评论中提到的，您本质上是截断或泄漏。

Answer 2

您可以像您一样使用傅里叶位移定理来做到这一点，但需要添加过滤。我问了类似的问题here。结果看起来“错误”的原因是因为您的输入向量不连续。您得到的结果确实是“正确的”（或至少是true）。

您看到的问题称为“吉布斯振铃”。您可以通过使用阶跃响应中没有振铃的低通滤波器（this维基百科链接很好地解释了该解决方案）来减轻这种极端情况。尝试使用和不使用高斯滤波器的代码。这种伪影经常出现在 MRI 成像（以及许多其他信号处理情况）中，有时可以通过滤波来缓解。

Answer 3

import matplotlib.pyplot as plt
import numpy as np
 
def plotRealImag(x, y, xnew, ynew, title):
    plt.figure()
    plt.plot(x, np.real(y), 'o-', x, np.imag(y), 'o-')
    plt.plot(xnew, np.real(ynew), '.-', xnew, np.imag(ynew), '.-')
    plt.legend(['y real', 'y imag', 'y new real', 'y new imag'], loc='best')
    plt.title(title)
     
def plotMagPhase(x, y, xnew, ynew, title):
    plt.figure()
    plt.plot(x, np.abs(y), 'o-')
    plt.plot(xnew, np.abs(ynew), '.-')
    plt.legend(['y', 'y new'], loc='best')
    plt.title(title + " Magnitude")
     
    plt.figure()
    plt.plot(x, np.angle(y), 'o-')
    plt.plot(xnew, np.angle(ynew), '.-')
    plt.legend(['y', 'y new'], loc='best')
    plt.title(title + " Phase")
 
def sampleShift(y, shift):
    n = y.shape[0]
    Y = np.fft.fft(y)
    w = np.arange(n)/n
    w[-n//2:] -= 1 # Make phase ramp continuous across DC
    Ynew = np.exp(-2j*np.pi*w*shift)*Y
    ynew = np.fft.ifft(Ynew)
    return ynew
     
# Create data
n = 16
u = 10
x = np.linspace(0, 10, n, endpoint=False)
y = np.cos(-x**2/6.0) + 1j*np.sin(-x**2/6.0)
xnew = np.linspace(0, 10, u*n, endpoint=False)
 
# Shift data
shift = 0.3
ynew = sampleShift(y, shift)
delx = x[1]-x[0]
plotRealImag(x, y, x-shift*delx, ynew, "FFT Sample Shift")
plotMagPhase(x, y, x-shift*delx, ynew, "FFT Sample Shift")
 
plt.show()