set(0,'DefaultLineLineWidth',2);
set(0,'DefaultStemLineWidth',2);
figure(1);

fp = 0.2;
fs = 0.3;
fc = 0.5*(fp + fs);
df = fs - fp;
N = 15; % Type 1 LPF

% Check
if(ceil(fp*N) > floor(fs*N))
    fprintf('Change N\n');
    return;
end

% Helper
NFFT = 2^16;
freq = [0:NFFT-1]/NFFT;
ipass = find(freq<=fp);
istop = find(fs<=freq & freq<=0.5);
ibin = [ipass,istop];

% Design
L = N; % Length
M = L-1; % Order
g = M/2; % Group delay
k = [0:g]; % Frequency bins
f = k/N; % Frequencies [cycles/sample]
[~,minind]=sort(abs(f-fc),'ascend');
minind = minind(1);
A = double(f <= fc); % LPF magnitude response
num = 20;
vals = zeros(num,6);
ivals = 0;
for a = 0:1/num:1
    A(minind) = a;
    H = A .* exp(-1j*2*pi*g*f); % Add linear phase response
    Hs = [H, conj(H(end:-1:2))]; % Hermitian symmetry
    h = real(ifft(Hs));
    % Check design
    H = abs(fft(h,NFFT)); % Compute magnitude response
    [Ep,ipmax] = max(abs(1-H(ipass))); % Maximum pass band error
    [Es,ismax] = max(H(istop)); % Maximum error stop band error
    ivals = ivals + 1;
    vals(ivals,:) = [minind, a, Ep, freq(ibin(ipmax)), Es, freq(ibin(ismax))];
    disp(vals(ivals,:));
    subplot(311);
    plot(freq,H); hold on;
    ax = axis;
    plot(fp*[1 1],ax(3:4),'r');
    plot(fs*[1 1],ax(3:4),'r');
    stem(f,abs(Hs(1:g+1)));
    hold off;
    xlim([0 0.5]);
    subplot(312);
    plot(freq,20*log10(H)); hold on;
    ax = axis;
    plot(fp*[1 1],ax(3:4),'r');
    plot(fs*[1 1],ax(3:4),'r');
    hold off;
    xlim([0 0.5]);
    ylim([-70 10]);
    subplot(313);
    stem([0:M],h);
    shg;
    pause;
end
figure(2);
plot(vals(:,2),vals(:,[3,5]));


