fp = 0.2;
fs = 0.3;
wp = fp*2*pi;
ws = fs*2*pi;
wc = 0.5*(wp+ws);

L = 41;
M = L-1;
g = M/2;

b = zeros(g+1,1);
A = zeros(g+1,g+1);
D = @(w) double(w <= wc);
W = @(w) double(w <= wc) + 100*double(w > wc);
for k=0:g
    I = @(w) W(w) .* D(w) .* cos(w*k);
    b(k+1) = quadgk(I,0,wp) + quadgk(I,ws,pi);
    for m=0:g
        I = @(w) W(w) .* cos(w*k).*cos(w*m);
        A(k+1,m+1) = quadgk(I,0,wp) + quadgk(I,ws,pi);
    end
end

a = inv(A)*b;
h = [0.5*a(end:-1:2); a(1); 0.5*a(2:end)];

set(0,'DefaultLineLineWidth',2);
NFFT = 2^14;
f = [0:NFFT-1]/NFFT - 0.5;
plot(f,20*log10(abs(fftshift(fft(h,NFFT)))));
grid on;
shg;
