% FIR filter Design
  1. Frequency Sampling
  2. Window
  3. Optimal (Parks-McClellan / Remez)


Windowing method (fir1, fir2, kaiserord)

1) Design H(f): ideal low-pass filter
2) Inverse DFT(H(f)) --> h(t)
3) Evaluate h(t)

              M
W(exp(iw)) = sum  exp(-iwn) = (1 - exp(-iw(M+1)) / (1 - exp(-iw))
             n=0

           = exp(-iw(M/2) ((sin(w(M+1)/2) / sin(w/2))

[H,W,s]=freqz(1);
s.xunits = '';
s.yunits = 'dB';
s.plot='mag';

w = boxcar(32);  %rectangle
y = fft(w, 256);
x = real(y);
x1 = [x(128:256); x(1:127)];
plot(x1);

t = [0:255];

%Truncated Sinc Lowpass FIR filter
b = 0.4 * sinc(0.4 * (-25:25));
subplot(311), plot(b);
grid;
axis tight;
[H, W] = freqz(b, 1, 1024, 2);
s.yunits='db';
subplot(312), freqzplot(H, W, s);
s.yunits='squared';
subplot(313), freqzplot(H, W, s);

% Rectangle window, different number of taps
figure;
s.yunits='squared';
b = 0.4 * sinc(0.4 * (-25:25));
[H1, W] = freqz(b, 1, 1024, 2);
subplot(311), freqzplot(H1, W, s);
axis([0.2, 0.6, 0, 1.4]);

b = 0.4 * sinc(0.4 * (-50:50));
[H2, W] = freqz(b, 1, 1024, 2);
subplot(312), freqzplot(H2, W, s);
axis([0, 1, -100, 20]);
axis([0.2, 0.6, 0, 1.4]);

b = 0.4 * sinc(0.4 * (-200:200));
[H3, W] = freqz(b, 1, 1024, 2);
subplot(313), freqzplot(H3, W, s);
axis([0.2, 0.6, 0, 1.4]);

% 101-taps different windows
figure;
s.yunits='squared';
b = 0.4 * sinc(0.4 * (-50:50));
[H1, W] = freqz(b, 1, 1024, 2);
subplot(311), freqzplot(H1, W, s);
axis([0.2, 0.6, 0, 1.4]);

b = 0.4 * sinc(0.4 * (-50:50));
b = b .* bartlett(101)';
[H2, W] = freqz(b, 1, 1024, 2);
subplot(312), freqzplot(H2, W, s);
axis([0, 1, -100, 20]);
axis([0.2, 0.6, 0, 1.4]);

b = 0.4 * sinc(0.4 * (-50:50));
b = b .* hamming(101)';
[H3, W] = freqz(b, 1, 1024, 2);
subplot(313), freqzplot(H3, W, s);
axis([0.2, 0.6, 0, 1.4]);

% 101-taps more different windows
figure;
s.yunits='squared';
b = 0.4 * sinc(0.4 * (-50:50));
b = b .* hann(101)';
[H1, W] = freqz(b, 1, 1024, 2);
subplot(311), freqzplot(H1, W, s);
axis([0.2, 0.6, 0, 1.4]);

b = 0.4 * sinc(0.4 * (-50:50));
b = b .* blackman(101)';
[H2, W] = freqz(b, 1, 1024, 2);
subplot(312), freqzplot(H2, W, s);
axis([0, 1, -100, 20]);
axis([0.2, 0.6, 0, 1.4]);

b = 0.4 * sinc(0.4 * (-50:50));
b = b .* hamming(101)';
[H3, W] = freqz(b, 1, 1024, 2);
subplot(313), freqzplot(H3, W, s);
axis([0.2, 0.6, 0, 1.4]);


% 151-taps different window Magnitude plots
figure;
s.yunits = 'dB';
% 151-tap FIR different windows
b = fir1(151,[0.1], boxcar(152)); %lowpass
[H,W] = freqz(b,1,512);
h = subplot(411);;
freqzplot(H, W, s);
axis([0, 0.6, -80, 10]);

b = fir1(151,[0.1], bartlett(152)); %lowpass
[H,W] = freqz(b,1,512);
h = subplot(412);;
freqzplot(H, W, s);
axis([0, 0.6, -80, 10]);

b = fir1(151,[0.1]); %lowpass  Hamming
[H,W] = freqz(b,1,512);
h = subplot(413);;
freqzplot(H, W, s);
axis([0, 0.6, -80, 10]);

b = fir1(151,[0.1], kaiser(152, 4)); %lowpass
[H,W] = freqz(b,1,512);
h = subplot(414);;
freqzplot(H, W, s);
axis([0, 0.6, -80, 10]);

% fir1 hamming different number of taps
s.yunits='dB';
b1 = fir1(5,[0.1]); %lowpass
[H,W] = freqz(b1,1,512);
h = subplot(411);;
freqzplot(H, W, s);
axis([0, 2, -80, 10]);

b2 = fir1(15,[0.1]); %lowpass
[H,W] = freqz(b2,1,512);
subplot(412), freqzplot(H, W, s);
axis([0, 2, -80, 10]);

b3 = fir1(51,[0.1]); %lowhpass
[H,W] = freqz(b3,1,512);
subplot(413), freqzplot(H, W, s);
title('51-tap');
axis([0, 2, -80, 10]);

b4 = fir1(201,[0.1]); %lowpass
[H,W] = freqz(b4,1,512);
subplot(414), freqzplot(H, W, s);
title('201-tap');
axis([0, 2, -80, 10]);

figure, zplane(b1);
figure, zplane(b2);
figure, zplane(b3);
figure, zplane(b4);

% different types
b = fir1(48,[0.65], 'high'); %highpass
[H,W] = freqz(b,1,512);
subplot(311), freqzplot(H, W, s);
axis([0, 3.14, -100, 10]);

b = fir1(48,[0.35 0.65]); %bandnpass
[H,W] = freqz(b,1,512);
subplot(312), freqzplot(H, W, s);
axis([0, 3.14, -100, 10]);

b = fir1(48,[0.35 0.65],'stop'); %bandstop
[H,W] = freqz(b,1,512);
subplot(313), freqzplot(H, W, s);
axis([0, 3.14, -100, 10]);

% Multiband: fir2
n = 50;
f = [0 .2 .5 .6 .7 .8 1];
m = [1 0 0 1 1 0 0];
b = fir2(n,f,m);
freqz(b,1,512);

Optimal (firls, remez)
Least-Square
Parks-McClellan

h = remez(51, [0 0.125 0.25, 1], [ 1 1 0 0]);
H = fft(h, 1024);
mag = abs([H(513:1024), H(1:512)]);
plot([-pi:pi/512:pi-pi/1024], mag);
figure; freqz(h);

% specify ripple deviation
Fs=44100;
F=[3300 4000];
M=[1 0];
dev=[0.002 0.002];  % ripple deviation
[N Fo Mo W]=remezord(F, M, dev, Fs)
h=remez(N, Fo, Mo, W);
[H, W]=freqz(h, 1, 1024, Fs);
freqzplot(H, W, s);
axis([0, 8000, -100, 10]);


% Comparison of FIR design methods
N=100;
h=remez(N, [0 0.1 0.2 1], [1 1 0 0]);
[H, W]=freqz(h);
subplot(311), freqzplot(H, W, s);
axis([0, 1, -100, 10]);
title('REMEZ');

h=firls(N, [0 0.1 0.2 1], [1 1 0 0]);
[H, W]=freqz(h);
subplot(312), freqzplot(H, W, s);
axis([0, 1, -100, 10]);
title('Leaset-Square');

h=fir2(N, [0 0.1 0.2 1], [1 1 0 0]);
[H, W]=freqz(h);
subplot(313), freqzplot(H, W, s);
axis([0, 1, -100, 10]);
title('Windowing');

