离散时间序列的内插算法（利用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
上面的代码只是示意性的，没有考虑计算效率一类的问题。实际应用的话，可以在这个基础上进一步优化。
郑重声明：本站内容如果来自互联网及其他传播媒体，其版权均属原媒体及文章作者所有。转载目的在于传递更多信息及用于网络分享，并不代表本站赞同其观点和对其真实性负责，也不构成任何其他建议。
