% Chapter 11, Example 11.2, Fig.11.6
% Eric Dubois, updated 2019-01-14
% Plots for figure which shows three
% filters for 16:1 separable upsampling: bilinear, bicubic and a designed
% least-pth filter with frequency domain constraints
% Required functions: ideal_response_qs, filter_design_qs_con

clear all, close all
%interpolation factor is K in each dimension
K=4;
% Bilinear interpolator
n1 = [-(K-1):(K-1)];
h1 = 1 - abs(n1)/K;
h11 = h1'*h1;
% Perspective plot of frequency response
[H11,fx,fy] = freqz2(h11,64,64,[1 1]);
[Fx,Fy] = meshgrid(fx,fy);
%perspective plot
figure;

mesh(Fx(1:2:64,1:2:64),Fy(1:2:64,1:2:64),abs(H11(1:2:64,1:2:64)));            % generate the perspective plot
colormap([0 0 0]);     % use black only
xlabel('u (c/opx)'), ylabel('v (c/opx)');
set(gca,'ydir','reverse'); 
%write the plot to a file
set(gcf,'Color',[1 1 1]);

% Bicubic interpolator
rng1 = (0:K-1)/K;
bcp1 = 1.5*(rng1).^3 -2.5*(rng1).^2 +1;
rng2 = (K:2*K-1)/K;
bcp2 = -0.5*(rng2).^3 + 2.5*(rng2).^2 - 4*(rng2) + 2;
h3 = [fliplr(bcp2) fliplr(bcp1) bcp1(2:size(bcp1,2)) bcp2]; 
h33 = h3'*h3;
% Perspective plot of frequency response
[H33,fx,fy] = freqz2(h33,64,64,[1 1]);
[Fx,Fy] = meshgrid(fx,fy);
%perspective plot
figure;

mesh(Fx(1:2:64,1:2:64),Fy(1:2:64,1:2:64),abs(H33(1:2:64,1:2:64)));            % generate the perspective plot
colormap([0 0 0]);     % use black only
xlabel('u (c/opx)'), ylabel('v (c/opx)');
set(gca,'ydir','reverse'); 
%write the plot to a file
set(gcf,'Color',[1 1 1]);

%design an interpolation filter with quadrantal symmetry and constraints on
%the frequency response

%filter size is L x L
L = 29;
%generate the ideal response
% the ideal passband edge is 1/2K. The the passband edge be a1/2K and the
% stopband edge be a2/2K where a1 < 1 and a2 > 1.
a1 = 0.8;
a2 = 1.2;
pbe = a1/(2*K);
sbe = a2/(2*K);
% the ideal passband gain is G = K^2
G = K^2;
%generate a matrix of ideal response values at the origin, the passband and
%stopband edges and the edge of the unit cell.
spec(1,:) = [0 0 G];  % origin
%generate 5 values on each of the horizontal and vertical passband edges 
for i = 1:5
    spec(1+i,:) = [pbe,-pbe + (i-1)*pbe/2 G];
    spec(6+i,:) = [-pbe + (i-1)*pbe/2 pbe G];
end
%generate 5 values on each of the horizontal and vertical stopband edges 
for i = 1:5
    spec(11+i,:) = [sbe,-sbe + (i-1)*sbe/2 0];
    spec(16+i,:) = [-sbe + (i-1)*sbe/2 sbe 0];
end
%generate 5 values on each of the horizontal and vertical unit cell edges
for i = 1:5
    spec(21+i,:) = [0.5,-0.5 + (i-1)/4 0];
    spec(26+i,:) = [-0.5 + (i-1)/4 0.5 0];
end
%save the specs to temporary text file
dlmwrite('interp_qs.txt',spec,'delimiter',' ');
%generate a matrix of frequency domain constraints
%There are  K^2 constraints, but with quadrantal symmetry, there are
%(fix(K/2)+1)^2 independent constraints for K even
ind = 1;
for j=0:fix(K/2)
    for i=0:fix(K/2)
        if (i == 0) && (j == 0)
            const(ind,:) = [0 0 K^2];
            ind = ind+1;
        else
            
            const(ind,:) = [i/K j/K 0];
            ind = ind+1;
        end
    end
end
%save constraints to temporary text file
dlmwrite('interp_qs_cons.txt',const,'delimiter',' ');
%Design the interpolation filter
%generate the ideal response
%Generate ideal filter response
Npt=64;
[f1,f2] = freqspace(Npt);
[F1,F2] = meshgrid(.5*f1,.5*f2);
Hd = ideal_response_qs('interp_qs.txt',Npt);
%perspective plot of the ideal response
figure;
mesh(Fx(1:2:64,1:2:64),Fy(1:2:64,1:2:64),abs(Hd(1:2:64,1:2:64)));            % generate the perspective plot
%view(-18,40);
colormap([0 0 0]);     % use black only
xlabel('u (c/opx)'), ylabel('v (c/opx)');
set(gca,'ydir','reverse'); 
%write the plot to a file
set(gcf,'Color',[1 1 1]);
%Design filter and save the coefficients to a file
h_int = filter_design_qs_con(Hd,L,L,'interp_qs_cons.txt');
%save('data/interp16_qs_coeff.mat','h_int');
%Compute and plot the frequency response
[H,fx,fy] = freqz2(h_int,64,64,[1 1]);
[Fx,Fy] = meshgrid(fx,fy);
%perspective plot
figure;
mesh(Fx(1:2:64,1:2:64),Fy(1:2:64,1:2:64),abs(H(1:2:64,1:2:64)));            % generate the perspective plot
%view(-18,40);
colormap([0 0 0]);     % use black only
xlabel('u (c/opx)'), ylabel('v (c/opx)');
set(gca,'ydir','reverse'); 

set(gcf,'Color',[1 1 1]);