interpolation_by_dfs.m

转化工具

function [xi,Xi] = interpolation_by_DFS(T,x,Ws,ti)
%T : sampling interval (sample period)
%x : discrete-time sequence
%Ws: normalized stop frequency (1.0=pi[rad])
%ti: interpolation time range or # of divisions for T
if nargin < 4, ti = 5; end
if nargin < 3 | Ws > 1, Ws = 1; end
N = length(x);
if length(ti) == 1
ti = 0:T/ti:(N-1)*T; %subinterval divided by ti
end
ks = ceil(Ws*N/2);
Xi = fft(x);
Xi(ks + 2:N - ks) = zeros(1,N - 2*ks - 1); %filtered spectrum
xi = zeros(1,length(ti));
for k = 2:N/2
xi = xi+Xi(k)*exp(j*2*pi*(k - 1)*ti/N/T);
end
xi = real(2*xi+Xi(1)+Xi(N/2+1)*cos(pi*ti/T))/N; %Eq.(3.9.5)
%interpolate_by_DFS
clear, clf
w1 = pi; w2 = .5*pi; %two tones
N = 32; n = [0:N - 1]; T = 0.1; t = n*T;
x = sin(w1*t)+0.5*sin(w2*t)+(rand(1,N) - 0.5); %0.2*sin(20*t);
ti = [0:T/5:(N - 1)*T];
subplot(411), plot(t,x,'k.') %original data sequence
title('original sequence and interpolated signal')
[xi,Xi] = interpolation_by_DFS(T,x,1,ti);
hold on, plot(ti,xi,'r') %reconstructed signal
k = [0:N - 1];
subplot(412), stem(k,abs(Xi),'k.') %original spectrum
title('original spectrum')
[xi,Xi] = interpolation_by_DFS(T,x,1/2,ti);
subplot(413), stem(k,abs(Xi),'r.') %filtered spectrum
title('filtered spectrum')
subplot(414), plot(t,x,'k.', ti,xi,'r') %filtered signal
title('filtered/smoothed signal')