% M07-filt.diary
% Design a 3-pt filter to remove one of two sine tones mixed together
% Dan Ellis dpwe@ee.columbia.edu 
% 2007-10-11

% filter is [alpha beta alpha]
% analysis gives mag. resp. as beta + 2 alpha cos(w)
% so to make it zero at w = 0.4, we have
% alpha = -beta/(2 cos(0.4))
% and to make it one at w = 0.1, we have
% beta (1 + 2(-1/(2 cos(0.4)))cos(0.1)) = 1
% hence...
beta  = 1/(1 + 2*(-1/(2*cos(0.4)))*cos(0.1))
% beta =
%   -12.4563
alpha = -beta/(2*cos(0.4))
% al =
%     6.7619

% Plot results...
% First, the impulse response
h = [alpha beta alpha];
subplot(311)
stem(h)
axis([0 50 -15 15])
grid

% Then the magnitude response (DTFT)
subplot(312)
[H,W] = freqz(h);
% Plot the magnitude in dB (log scaling)
plot(W, 20*log10(abs(H)))
axis([0 1.2 -20 20]);
grid
% gain of zero at w = 0.4 corresponds -Inf dB
% notice that gain goes through 0 dB (=1x) just around w = 0.1

% Finally, compare signal before and after filtering
nn = 1:8000;
w1 = 0.4;
w2 = 0.1;
x = 0.5*(sin(w1*nn) + sin(w2*nn));
y = conv(h,x);
subplot(313)
plot(1:200,x(1:200),1:200,y(1:200),'r')
axis([0 200 -1 1])
grid
% Check out the sound
sound(x,16000);
sound(y,16000);
% Higher component (w1) clearly removed