interpolation in the time domain by zero Padding In The Frequency Domain

xinglong-liu

最近学习DSP时，在网上看到这样一个网页，作者也很出名，相信很多人都读过他的书, understanding digital signal processing
webpage title:

How to Interpolate in the Time-Domain by
Zero-Padding in the Frequency Domain-by Rick Lyons

http://www.dspguru.com/dsp/howtos/how-to-interpolate-in-time-domain-by-zero-padding-in-frequency-domain

根据他描述的方法，我用matlab编程验证了一下，发现结果有些出入，不知问题出在哪里。因此贴在这里希望大家能够给予指点。






function test()
close all
clear all
clc;
dt=1/8e3;
df=1/dt;
N=7;
t=[0:1:N-1]*dt;
f1=1e3;
f2=2e3;
s(:,1)=sin(2*pi*f1*t)+0.5*sin(2*pi*f2*t+0.75*pi);
% s=s.*hanning(len);
S=fft(s);
f=[0:1:N-1]*df/N;
figure(1)
subplot(211)
plot(t,s,'-O');hold on
figure(2)
subplot(211)
plot(real(S),'Ob');hold on
subplot(212)
plot(imag(S),'sr');hold on
M=2*N;
if mod(N,2)==0
    Sm=[S(1:N/2);zeros(M/2,1); S(1+N/2); zeros(M/2,1);S(2+N/2:N)];
else
    Sm=[S(1:(N+1)/2);zeros(M,1);S((N+3)/2:N)];
end
figure(3)
subplot(211)
plot(real(Sm),'Ob');hold on
subplot(212)
plot(imag(Sm),'sr');hold on
sn=ifft(Sm);
Nn=length(sn);
Re_sn=real(sn)*(N+M)/N;
Im_sn=imag(sn)*(N+M)/N;
dt=max(t)/(M+N-1);
tn=[0:M+N-1]*dt;
figure(1)
subplot(211)
plot(tn,Re_sn,'-*r');hold on
grid on
legend('orginal data','interpolated data')

[ 本帖最后由 xinglong-liu 于 2010-1-21 15:41 编辑 ]

songzy41

程序基本上是对的，只是计算内插后的dt和tn不正确，把
dt=max(t)/(M+N-1);
tn=[0:M+N-1]*dt;
改为
dt=dt*N/(N+M);
tn=(0:N+M-1)*dt;
就能得到正确的结果了。

xinglong-liu

谢谢songzy41的指点，终于得到了正确的结果！
还需要继续学习！