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

fp = 0.2;
fs = 0.3;
fc = 0.5*(fp + fs);
df = fs - fp;
N = 27; % 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');
mininda = minind(1);
minindb = minind(2);
A = double(f <= fc); % LPF magnitude response
num = 21;
vec = 0:1/(num-1):1;
vals = zeros(num,num,2);
Esmin = 1;
amin = 0;
bmin = 0;
for ia = 1:num
    for ib = 1:num
        A(mininda) = vec(ia);
        A(minindb) = vec(ib);
        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
        vals(ia,ib,1) = Ep;
        vals(ia,ib,2) = Es;
        disp([minind, a, Ep, freq(ibin(ipmax)), Es, freq(ibin(ismax))]);
        
        if 0
            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;
            %drawnow;
        end
        
        if(Es < Esmin)
            Esmin = Es;
            amin = vec(ia);
            bmin = vec(ib);
            %pause;
        end
    end
end

clf;
% Take the best one found during the search
A(mininda) = amin;
A(minindb) = bmin;
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

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;

% Add hamming windowed FIR
n = [0:M].';
h = (2*fc)*sinc(2*fc*(n-g)).*hamming(L);
H = abs(fft(h,NFFT));
subplot(311); hold on;
plot(freq,H); hold off;
subplot(312); hold on;
plot(freq,20*log10(H)); hold off;
subplot(313); hold on;
stem([0:M],h); hold off;
shg;



