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

利用FFT实现信号内插

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

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

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

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

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

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

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

当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 

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