%  dip07_2.m               gr23mar15   26oct17   nov20
%   interpolation (1-D)
clear all;

f = [3 5 7 7 2 5 2 3 6];    % data
return;

%%  --------- PSEUDO-CONTINUOUS interpolation using sinc, R0, R1
%%% see also dsp: Ex_2_2.m
%
figure(1);  stem([-4:4],f,'k','filled');  axis([-4,4,-2,8.5]);
hold on;  pause();
%				 passes through data values; 
x = -4:.01:4;   %   "continuous" interpolation: sample rate = .01
y = 0;
for i=-4:4;
    s = f(i+5) * sinc(x-i);   %  sinc, centered on i-th position
%    s = f(i+5) * double( (x-i) >= -.5 & (x-i) < .5 );  % R0(x-i);
%    s = f(i+5) * max( 1 - (x-i) .* sign(x-i), 0);  % R1(x-i);
    plot(x,s); pause(.5);
    y = y + s;
    end;
pause();
plot(x,y,'r','linewidth',3); hold off; 

%%  ---------- DISCRETE interpolation via upsampling and filtering
%	approximation (curve fitting) rather than interpolation,  if order >=1
ups = 7;   % integer upsampling factor, e.g. 3, 10, whatever
ord = 0;    % order of the interpolator (0,...3)
z = zeros(ups-1,length(f));
fu = [f; z];  fu = fu(:)';      % upsampled data
fu = fu(1:end-ups+1);           % truncate tail of zeros
h0 = ones(1,ups);       % zero-order interpolator
h = h0;
for i=1:ord             % build (ord)-order interpolator
    h = conv(h,h0/ups);
    end;
figure(2); stem(h);     % impulse response of interpolator
g = conv(fu,h);         % filter upsampled data
m = floor(length(h)/2);  g = g(m:end-m+1);      % truncate tails
figure(3);
stem(fu,'filled'); axis([1,length(g),-1,8.5]); hold on;
stem(g,'r');    hold off;

%%  ---------------------- Keys
step = 1 / ups;
x = -2+step : step : -1-step;  K1 = -.5*abs(x.^3)+2.5*x.^2-4*abs(x)+2;
x = -1 : step : 1;             K2 = 1.5*abs(x.^3)-2.5*x.^2+1;
x = 1+step : step : 2;         K3 = -.5*abs(x.^3)+2.5*x.^2-4*abs(x)+2;
h = [K1,  K2,  K3];
figure(2); stem(h);   % impulse response of interpolator
g = conv(fu,h);       % filter upsampled data
m = floor(length(h)/2);  g = g(m:end-m+1);      % truncate tails
figure(3);
stem(fu,'filled'); axis([1,length(g),-1,8.5]); hold on;
stem(g,'r');    hold off;

%%  ---------------------- R1  same as above when ord=1
step = 1 / ups;
x = -1+step : step : 0;  R11 = x + 1;
x =  0+step : step : 1;  R12 = 1 - x;
R1 = [R11, R12];
figure(2); stem(R1);   % impulse response of interpolator
g = conv(fu,R1);       % filter upsampled data
m = floor(length(R1)/2);  g = g(m:end-m+1);      % truncate tails
figure(3);
stem(fu,'filled'); axis([1,length(g),-1,8.5]); hold on;
stem(g,'r');    hold off;

%%                     Response in the DFT domain
%  compare to sinc interpolation
x = 0:10; 
y = sin(x) + cos(3*x);
figure(1); stem(x,y); hold on; pause();
xi = 0:.1:10; 
 yi = interp1(x,y,xi,'nearest');
% yi = interp1(x,y,xi,'linear');
% yi = interp1(x,y,xi,'cubic');

ysinc = 0;     % sinc interpol.
for i=0:10;    ysinc = ysinc + y(i+1)*sinc(xi-i);    end;
plot(xi,yi,xi,ysinc);
hold off

figure(2);
Yi = abs(fft(yi,1024)); 
Ysinc = abs(fft(ysinc,1024)); 
plot([0:512],Yi(1:513),[0:512],Ysinc(1:513)); axis('tight');
return;

%%  ---------- interp. by factor L/M (slide 14)
L = 5; M = 2;
order = 20;
h_up = L * fir1(order, 1/L);	% design LP filter, ome_c=pi/L, gain L
h_down = fir1(order, 1/M);	    % design LP filter, ome_c=pi/M, gain 1

x = filter(.1, [1 -.9], randn(1,1000)); % input signal: colored noise
xu = [x; zeros(L-1,length(x))];
xu = xu(:)';               % upsampled data
w = filter(h_up,1,xu);     % apply anti-imaging filter
ww = filter(h_down,1,w);   % apply anti-aliasing filter (unnecessary)
y = ww(1:M:end);           % downsample
 
figure(21); subplot(221); plot(abs(fft(x)));
            subplot(222); plot(abs(fft(xu)));
            subplot(223); plot(abs(fft(w)));
            subplot(224); plot(abs(fft(y)));
figure(25); subplot(211); plot(x);
            subplot(212); plot(y);


