% Chapter 11, Fig.11.10
% Eric Dubois, updated 2019-01-14
%Design a 2D filter for 2:1 downsampling in each direction (4:1 overall)
%Will apply it the Barbara image for Fig.11.10
% Required functions: ideal_response_qs, filter_design_qs_con
% Design a filter with quadrantal symmetry and DC gain of 1.0
clear all, close all
K=2; % downsampling factor in each dimension
%filter size is L x L
L = 15;
%generate the ideal response
% the ideal passband edge is 1/2K. Let 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 = 1
G = 1;
%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 1];
    spec(6+i,:) = [-pbe + (i-1)*pbe/2 pbe 1];
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('antialias_qs.txt',spec,'delimiter',' ');

%A single frequency domain constraint at DC
%save constraints to temporary text file
const(1,:) = [ 0 0 1];
dlmwrite('antialias_qs_cons.txt',const,'delimiter',' ');

%Design the interpolation filter
%Generate ideal filter response
Npt=64;
[f1,f2] = freqspace(Npt);
[F1,F2] = meshgrid(.5*f1,.5*f2);
Hd = ideal_response_qs('antialias_qs.txt',Npt);
%Design filter
h_aa = filter_design_qs_con(Hd,L,L,'antialias_qs_cons.txt');
save('antialias_qs_coeff.mat','h_aa');
%Compute and plot the frequency response
[H_aa,fx,fy] = freqz2(h_aa,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_aa(1:2:64,1:2:64)));            % generate the perspective plot
%view(-18,40);
colormap([0 0 0]);     % use black only
%axis([-.5 .5 -.5 .5 0 1.2]);
xlabel('u (c/ipx)'), ylabel('v (c/ipx)');
set(gca,'ydir','reverse');