% Chapter 10, Fig. 10.13
%Eric Dubois, updated 2018-12-20
%Circularly symmetric ideal response examples
%One-dimensional radial profile of several ideal responses with a 3dB
%attentuation at 0.125 c/px
clear all, close all;
three_db = 0.125;
ur = 0:.001:.5; % frequencies at which to sample the radial profile
%"ideal" low pass filter, step function in frequency domain
Hr_step = ones(size(ur));
Hr_step(ur>three_db) = 0;
%linear transition band
width = 0.15; % width of the transition band
Hr_lin = -ur/width + 1/sqrt(2) + three_db/width;
Hr_lin = min(1,Hr_lin);
Hr_lin = max(0,Hr_lin);
%
%raised-cosine transition band
a = three_db - acos(sqrt(2)-1)*width/pi;
Hr_raised_cos = 0.5*(1+ cos(pi*(ur - a)/width));
Hr_raised_cos(ur<a) = 1;
Hr_raised_cos(ur > a + width) = 0;
%Plot
H(1,:) = Hr_step;
H(2,:) = Hr_lin;
H(3,:) = Hr_raised_cos;
%Set plot size, default font and darker colors
figure

%set the axis limits prior to calling DashLine
axis([0 0.5 0 1]);
axdata = axis;

plot(ur,H(1,:),ur,H(2,:),ur,H(3,:));
%Labels
xlabel('Spatial frequency (c/px)');
ylabel('Ideal radial frequency response')
set(gcf,'Color',[1 1 1]);
set(gca,'fontname','times')          %set the plot font to Times

% Perspective plots of two-dimensional ideal response
%use 32 x 32 arrays to get a suitable density of lines
ux = -.5:1/32:.5; uy = ux;
[Ux,Uy] = meshgrid(ux,uy);
%Step transition
H_step = ones(size(Ux));
H_step(Ux.^2 + Uy.^2 > three_db^2) = 0;
figure;

mesh(Ux,Uy,H_step)
colormap([0 0 0]);     % use black only
xlabel('u (c/px)'), ylabel('v (c/px)');
set(gca,'ydir','reverse'); 
set(gca,'fontname','times')          %set the plot font to Times

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

%Linear transition
H_lin = -sqrt(Ux.^2 + Uy.^2)/width + 1/sqrt(2) + three_db/width;
H_lin = min(1,H_lin);
H_lin = max(0,H_lin);
figure;
mesh(Ux,Uy,H_lin);
colormap([0 0 0]);     % use black only
xlabel('u (c/px)'), ylabel('v (c/px)');
set(gca,'ydir','reverse'); 
set(gca,'fontname','times')          %set the plot font to Times
set(gcf,'Color',[1 1 1]);

%Raised cosine transition
Hr_raised_cos = 0.5*(1+ cos(pi*(sqrt(Ux.^2+Uy.^2) - a)/width));
Hr_raised_cos(sqrt(Ux.^2+Uy.^2)<a) = 1;
Hr_raised_cos(sqrt(Ux.^2+Uy.^2) > a + width) = 0;
figure;
mesh(Ux,Uy,Hr_raised_cos);
colormap([0 0 0]);     % use black only
xlabel('u (c/px)'), ylabel('v (c/px)');
set(gca,'ydir','reverse'); 
set(gca,'fontname','times')          %set the plot font to Times
set(gcf,'Color',[1 1 1]);