%% For plotting frequency response
NFFT = 2^14;
freq = [0:NFFT-1]/NFFT;

%% Type 1 LPF
M = 20; % Even order
L = M+1; % Odd length
nd = M/2; % Group delay
fc = 0.23; % Cut off frequency (normalized)

% Truncated ideal filter
n = [0:M].';
hd = 2*fc*sinc(2*fc*(n-nd));

% 1. Choose window shape to achieve
% stop band attenuation and passband
% ripple.
% 2. Choose filter order M to achieve
% transition band width.

win = my_rect(L);
%win = my_triang(L);
%win = my_hamming(L);
%win = my_hann(L);
%win = window(@rectwin        ,L);
%win = window(@bartlett       ,L);
%win = window(@triang         ,L);
%win = window(@barthannwin    ,L);
%win = window(@blackman       ,L);
%win = window(@blackmanharris ,L);
%win = window(@bohmanwin      ,L);
%win = window(@chebwin        ,L);
%win = window(@gausswin       ,L);
%win = window(@hamming        ,L);
%win = window(@hann           ,L);
%win = window(@kaiser         ,L);

% Apply window
h = hd.*win;

% Check frequency response
plot(freq,20*log10(abs(fft([hd h],NFFT))));
grid on;

% Classroom computer can't do this
%% Use Matlab built-ins
wp = 0.25*pi; fp = wp/(2*pi);
ws = 0.35*pi; fs = ws/(2*pi);
Fs = 1; % sample rate
Ap = 0.1; Bp = 10^(Ap/20); dp = (Bp-1)/(Bp+1);
As = -50; ds = 10^(As/20); 

if 0
[M,Wn,beta,filtype]=kaiserord([fp fs],[1 0],[dp ds],Fs);
if(mod(M,2)==1) M=M+1; L=M+1; end % M is odd, make it even
h = fir1(M,Wn,filtype,kaiser(L,beta),'noscale');
end

% Check frequency response
plot(freq,20*log10(abs(fft(h,NFFT))));
hold on;
ax = axis;
plot([0 fp],Ap*[1 1],'r');
plot([0 fp],-Ap*[1 1],'r');
plot([fs 0.5],As*[1 1],'r');
plot(fp*[1 1],ax(3:4),'r');
plot(fs*[1 1],ax(3:4),'r');
hold off;
grid on;
ylim([-100 10]);

if  0
%% Equiripple filter
[M,F,A,W] = firpmord([fp fs],[1 0],[dp ds],1);
h = firpm(M,F,A,W);
%h = firpm(M,[0 fp fs 0.5]/0.5,[1 1 0 0],[1 3]);
end
% Check frequency response
plot(freq,20*log10(abs(fft(h,NFFT))));
hold on;
ax = axis;
plot([0 fp],Ap*[1 1],'r');
plot([0 fp],-Ap*[1 1],'r');
plot([fs 0.5],As*[1 1],'r');
plot(fp*[1 1],ax(3:4),'r');
plot(fs*[1 1],ax(3:4),'r');
hold off;
grid on;
ylim([-100 10]);


%% Run signals through the filters
NFFT = 2^14; freq = [0:NFFT-1]/NFFT;
n = [0:1e4-1];
x = cos(2*pi*0.05*n) + ...
    cos(2*pi*0.15*n) + ...
    cos(2*pi*0.25*n) + ...
    cos(2*pi*0.26*n) + ...
    cos(2*pi*0.27*n);
y = conv(x,h,'valid');
X = abs(fft(x(:),NFFT));
Y = abs(fft(y(:),NFFT));
H = abs(fft(h(:),NFFT));
plot(freq,20*log10(X),'b'); hold on;
plot(freq,20*log10(Y),'r');
plot(freq,20*log10(H),'Color',[0 0.5 0],'LineWidth',2);
grid on;
hold off;
ylim([-100 80]);


    
    
