%FILTER_DESIGN_QS_UNC Design a 2-D filter with quadrantal symmetry by using
%   the minimax error between the ideal frequency response and the
%   approximation given by the window method as the design criteria. UNC
%   stands for unconstrained optimization.
%
%   h = filter_design_qs_unc(Hd,Nfilt_row,Nfilt_col)
%   Hd: the ideal frequency response
%   Nfilt_row: vertical order of the output filter
%   Nfilt_col: horizontal order of the output filter
%   h: output filter
%
%   This function will also print out the minimax error between the output
%   filter and the ideal filter.
%
%   Example
%   --------
%   Design a 15x15 lowpass filter based on the ideal frequency response Hd
%
%   h=filter_design_qs_unc(Hd,15,15);
%
%   NOTE: Hd can be designed using the function ideal_response_qs. Please
%   refer to the help documents of ideal_response_qs for more information.

function h = filter_design_unc(Hd,Nfilt_row,Nfilt_col)

Npt=size(Hd,1);
%Extract one quadrant of ideal filter response
Hd_qua = Hd((size(Hd,1)*0.5+1):size(Hd,1),(size(Hd,2)*0.5 + 1):size(Hd,2));
%Generate desired window
win=ones(Nfilt_row,Nfilt_col);
%Generate actual filter response using window method
ha = fwind2(Hd,win);
ha_qua0 = ha(((size(ha,1)-1)*0.5+1):size(ha,1) , ((size(ha,2)-1)*0.5+1):size(ha,2));
ha_vec0 = reshape(ha_qua0,((Nfilt_row-1)*.5+1)*((Nfilt_col-1)*.5+1),1);
%Set optimizing options
options = optimset('MaxFunEvals',100000,'TolFun',1e-15,'TolCon',1e-15,'TolX',1e-15);

%Optimize solution with constraints
%Find minimum-error solution
p=2;
exitcond=0;
sol = fminunc(@(ha_vec) p_error(ha_vec,Hd_qua,Npt,Nfilt_row,Nfilt_col,2),ha_vec0);
error_minimax=minmax_error(sol,Nfilt_row,Nfilt_col,Hd_qua,Npt);
while not(exitcond)
    p=p+2;
    sol_tem = fminunc(@(ha_vec) p_error(ha_vec,Hd_qua   ,Npt,Nfilt_row,Nfilt_col,p),ha_vec0);
    err_tem=minmax_error(sol_tem,Nfilt_row,Nfilt_col,Hd_qua,Npt);
    if (err_tem>=error_minimax)
        exitcond=1;
        p=p-2
    else
        sol=sol_tem;
        error_minimax=err_tem;
    end
end

%Display minimax error
error_minimax
%Convert solution vector to a quadrant
solqua=reshape(sol,(Nfilt_row-1)*0.5+1,(Nfilt_col-1)*0.5+1);
%Convert quadrant to full response
bottom = [fliplr(solqua(:,2:size(solqua,2))) solqua];
h = [flipud(bottom(2:size(bottom,1),:)) ; bottom];


%%%%%%%%%%%%%SUBFUNCTIONS to be used in this design%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Subfunction1: minmax_error
function err = minmax_error(vec,Nfilt_row,Nfilt_col,Hd_qua,Npt)
%Generate full impulse response and frequency response from one quadrant.
ha_qua=reshape(vec,(Nfilt_row-1)*0.5+1,(Nfilt_col-1)*0.5+1);
r=size(ha_qua,1);
c=size(ha_qua,2);
upright = flipud(ha_qua(2:r,:));
upleft = fliplr(upright(:,2:c));
dleft = fliplr(ha_qua(:,2:c));
ha_pad = [ha_qua zeros(r,Npt-2*c+1) dleft;zeros(Npt-2*r+1,Npt);
    upright zeros(r-1,Npt-2*c+1) upleft];
tem = fft2(ha_pad);
Ha_qua = real(tem(1:Npt/2,1:Npt/2));
%Calculate the minimax error
err = max(max(abs(Ha_qua - Hd_qua)));

%Subfunction2: p_error
%This function computes the p-error difference between the desired filter
%FREQUENCY response and the actual filter FREQUENCY response for one
%quadrant. 

%Input: ha_qua is the coefficient matrix for one quadrant of the actual 
%impulse response; Hd_qua is one quadrant of the ideal frequency response;
%Npt is the number of points on each dimension of Hd; Nfilt is the size of
%the actual impulse response; p is the power of the difference.
%Output:a scalar quantity indicating the pth power difference between 
%the actual frequency response and the ideal frequency response for one
%quadrant.
function err = p_error(ha_vec,Hd_qua,Npt,Nfilt_row,Nfilt_col,p)

%Generate full impulse response and frequency response from one quadrant.
ha_qua = reshape(ha_vec,(Nfilt_row-1)*0.5+1,(Nfilt_col-1)*0.5+1);
r=size(ha_qua,1);
c=size(ha_qua,2);
upright = flipud(ha_qua(2:r,:));
upleft = fliplr(upright(:,2:c));
dleft = fliplr(ha_qua(:,2:c));
ha_pad = [ha_qua zeros(r,Npt-2*c+1) dleft;zeros(Npt-2*r+1,Npt);
    upright zeros(r-1,Npt-2*c+1) upleft];
tem = fft2(ha_pad);
Ha_qua = real(tem(1:Npt/2,1:Npt/2));

%Calculate the p-error
err = sum(sum((Ha_qua - Hd_qua).^p));