%% Example 10.2 page 565

NFFT = 2^14;
freq = [0:NFFT-1]/NFFT - 0.5;

wp = 0.25*pi;
ws = 0.35*pi;
Ap = 0.1; % dB
As = 50; % dB

Bp = 10^(Ap/20);
dp = (Bp-1)/(Bp+1); % pass band ripple
Bs = (1+dp)*10^(-As/20);
ds = Bs; % stop band ripple (approx.)

%% Windowed filter design
Dw = ws - wp;
wc = 0.5*(ws + wp);
[d,dind] = min([dp ds]);
if(dind == 1)
    A = Ap; % design is dominated by pass band ripple
else
    A = As; % design is dominated by stop band ripple
end;
% Choose hamming window because it offers 53 dB attenuation in the
% stop band.
L = ceil(6.6*pi/Dw);
% Increase L to a type 1 filter
L = 67;
% Now we know the shape (Hamming) and the length.
win = hamming(L);
% Generate the impulse response
M = L-1; % filter order
alpha = M/2; % group delay
n = [0:L-1].';
hd = (wc/pi)*sinc((wc/pi)*(n-alpha));
h = hd.*win;
hhamming = h;

plot(freq,20*log10(abs(fftshift(fft(hd,NFFT)))),'Color',[0.8 0 0.3],'LineWidth',3);
hold on;
plot(freq,20*log10(abs(fftshift(fft(h ,NFFT)))),'Color',[0.0 0.8 0.3],'LineWidth',3);
ax = axis;
plot((wp/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
plot(-(wp/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
plot((ws/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
plot(-(ws/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
hold off;
grid on;
xlabel('Normalized Frequency [cyles/sample]');
ylabel('Magnitude [dB]');
legend('Non-Windowed','Windowed');
ylim([-100 10]);


%% Kaiser design
[M,Wn,beta,filtype]=kaiserord([wp ws]/pi/2,[1 0],[dp ds],1)
M = 60;
L = M+1;
alpha = M/2;
win = kaiser(L,beta);
n = [0:L-1].';
hd = (wc/pi)*sinc((wc/pi)*(n-alpha));
h = hd.*win;
hkaiser = h;

plot(freq,20*log10(abs(fftshift(fft(hd,NFFT)))),'Color',[0.8 0 0.3],'LineWidth',3);
hold on;
plot(freq,20*log10(abs(fftshift(fft(h ,NFFT)))),'Color',[0.0 0.8 0.3],'LineWidth',3);
ax = axis;
plot((wp/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
plot(-(wp/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
plot((ws/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
plot(-(ws/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
hold off;
grid on;
xlabel('Normalized Frequency [cyles/sample]');
ylabel('Magnitude [dB]');
legend('Non-Windowed','Windowed');
ylim([-100 10]);

%% Equiripple filter design
[M,F,A,W]=firpmord([wp ws]/pi/2,[1 0],[dp ds],1)
%M = M + 2;
h = firpm(M,F,A,W);
hfirpm = h;

plot(freq,20*log10(abs(fftshift(fft(hhamming,NFFT)))),'Color',[0.8 0 0.3],'LineWidth',3);
hold on;
plot(freq,20*log10(abs(fftshift(fft(hkaiser ,NFFT)))),'Color',[0.0 0.8 0.3],'LineWidth',3);
plot(freq,20*log10(abs(fftshift(fft(hfirpm  ,NFFT)))),'Color',[0.0 0.3 0.8],'LineWidth',3);
ax = axis;
plot((wp/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
plot(-(wp/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
plot((ws/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
plot(-(ws/(2*pi))*[1 1],ax(3:4),'Color',[0 0.5 0.5],'LineWidth',2);
hold off;
grid on;
xlabel('Normalized Frequency [cyles/sample]');
ylabel('Magnitude [dB]');
legend('Hamming','Kaiser','FIR-PM');
ylim([-100 10]);

fprintf('Hamming = %d\nKaiser  = %d\nFIR-PM  = %d\n',...
        length(hhamming),length(hkaiser),length(hfirpm));