离散时间序列的内插算法（利用fft）

有些时候，为了后续处理更方便，我们需要对采集到的数据点进行内插处理，也就是所谓的增采样。本文就来讨论一下常用的几种内插算法。

利用FFT实现信号内插

我们的信号 x(t) 是个实信号，带宽有限，能量有限。x
=x(nΔ)和 x’
=x(nΔ’)是对这个信号的两种采样，并且都满足采样定理的要求，也就是说信息并没有丢失。两次采样的采样率满足如下关系。



也就是说第二种采样的采样率是第一种采样的采样率的M 倍。如果我们有了x’
，那么很容易获得x
：



但是反过来就不是那么容易了。下面就来推导如何在已知x
的情况下，计算出x’
。

设x
序列的长度为 N，x’
序列的长度为 M×N。那么这两个序列的频谱分别如下：



同时，还有两个反变换公式：



从上面的推导我们首先看出，

，也就是通过这两个序列计算的频点是相同的。而这两个序列都满足乃奎斯特采样定律，因此对应的频谱上同频率的点的值应该是有关系的。而这个关系可以通过反变换公式来推导出。



可以看出，如果，则 x
= x’[nM]。但是这样求得的x’
会有许多元素是复数值。我们知道实数序列的傅立叶变换满足如下条件：

。为了扩展后还满足这个条件，可以这样做。



当N为奇数时：



当 N为偶数时：



这里要特别注意，N/2 那个元素需要分成两份，前后各放一半，只有这样变换后的结果才是对的。类似的处理在利用 FFT 计算离散解析信号时也会遇到。

下面给个scilab 的代码，scilab与 matlab 类似，数组下标从1开始，所以上面的的公式需要相应的修改一点。

function y = fft_interp(x, n)
    s = size(x);
    if s(1) == 1 then 
        nx = s(2);
        ny =  (nx) * n;
        s(2) = ny;
    else
        nx = s(1);
        ny =  (nx) * n;
        s(1) = ny;
    end    
    y = zeros(s(1), s(2));
    
    xx = fft(x);
    if (nx / 2) == int(nx / 2) then // 偶数
        s = nx/2;
        y(1:(1+s)) = xx(1:(1+s))
        y(ny-s+1:ny) = xx(nx-s+1:nx) 
        y(1+s) = y(1+s) / 2
        y(ny-s+1) = y(ny-s+1)/2
    else
        s = (nx - 1)/2;
        y(1:(1+s)) = xx(1:(1+s))
        y(ny-s+1:ny) = xx(s+2:nx)        
    end
    y = ifft(n*y)
endfunction

上面的代码只是示意性的，没有考虑计算效率一类的问题。实际应用的话，可以在这个基础上进一步优化。