% =====================================================
% ECE503 example problems for Lecture 2
% DRB 23-Jan-2012
%
% Note: Don't just run this Matlab code all at once.
% It is used to verify answers to multiple problems. 
% You should cut and paste pieces of this code into the 
% Matlab command window.
% =====================================================


% =====================================================
% EXAMPLE: Verify our solution to the zero-input and zero-state
% responses of the system defined by
% y[n] = 0.5y[n-1]+0.25y[n-2]+x[n]+0.5x[n-1]
% with y[-1] = 1 and y[-2] = 0
% and with input x[n] = 0.9^n u[n]
%
% You will probably need the class notes to understand
% where all the numbers came from for this example.
% =====================================================

% PRELIMINARIES
n = 0:100;      % sample indices
lam = roots([1 -0.5 -0.25]); % roots of characteristic polymomial
A = [1 1;lam'];              % form matrix to solve simultaneous equations

% Compute the zero-input, zero-state, and impulse responses using Matlab's 
% filter function 
yzim = filter([1 0.5],[1 -0.5 -0.25],zeros(1,length(n)),[0.5 0.25]); % zero-input
yzsm = filter([1 0.5],[1 -0.5 -0.25],0.9.^n,[0 0]); % zero-state
hm = filter([1 0.5],[1 -0.5 -0.25],[1 zeros(1,length(n)-1)],[0 0]); % impulse

% Compute our own solution to the zero-input response
% using the intitial conditions
alphazi = inv(A)*[0.5;0.5];     % solve simultaneous equations for y[0] and y[1]
yzi = alphazi(1)*lam(1).^n + alphazi(2)*lam(2).^n; % zero-input response

% Compute our own solution to the zero-state response
beta = (0.9^2+0.5*0.9)/(0.9^2-0.5*0.9-0.25);
alphazs = inv(A)*[1-beta; 1.9-0.9*beta]; % solve simultaneous equations
yzs = alphazs(1)*lam(1).^n + alphazs(2)*lam(2).^n+beta*0.9.^n; % zero-state response

% Compute our own solution to the impulse response
alphah = inv(A)*[1;1];      % solve simultaneous equations
h = alphah(1)*lam(1).^n + alphah(2)*lam(2).^n; % impulse response

% MAKE PLOTS
subplot(4,1,1)
stem(n,yzi);
hold on
stem(n,yzim,'r');
hold off
title('zero-input solution');
subplot(4,1,2)
stem(n,beta*0.9.^n); % particular solution
hold on
stem(n,yzsm,'r');
hold off
title('particular solution vs. zero-state response');
subplot(4,1,3)
stem(n,yzs);
hold on
stem(n,yzsm,'r');
hold off
title('zero-state response');
subplot(4,1,4);
stem(n,h);
hold on
stem(n,hm,'r');
hold off
title('impulse response');

% =====================================================
% EXAMPLE: Simple filtering. Find an FIR filter
% h = [alpha0,alpha1,alpha0] that blocks the tone
% at omega0 and passes the tone at omega1 with unit gain
% This follows the procedure in your textbook and 
% summarized in the slides.
% =====================================================
fs = 44100;                 % sampling frequency (Hz)
T = 5;                      % seconds of sound to generate
omega0 = (300/fs)*2*pi;
omega1 = (1000/fs)*2*pi;
n = 0:T*fs;
c0 = 1;
c1 = 1;
x = c0*cos(omega0*n)+c1*cos(omega1*n);
sound(x/2,fs);
disp('Press a key to filter out the low tone');
pause
alpha0 = 1/(2*(cos(omega1)-cos(omega0)));
alpha1 = 1-2*alpha0*cos(omega1);
h = [alpha0 alpha1 alpha0];
y = filter(h,1,x);
sound(y/2,fs);
[h,w]=freqz(h,1,1024);
plot(w,10*log10(abs(h).^2),[omega0 omega0],[-80 0],[omega1 omega1],[-80 0]);
xlabel('normalized freq (rad/sample)');
ylabel('magnitude response of the simple filter (dB)');
grid on

