% zeros on unit circle
z0 = -1;         % angle pi
z1 = 1;         % angle 0
numer = [1 -(z0+z1) z0*z1];
% poles inside unit circle
f0 = 440; fs = 8000;
w0hat = 2*pi*f0/fs;
r = 0.99;
p0 = r * exp(j*w0hat);
p1 = r * exp(-j*w0hat);
denom = [1 -(p0+p1) p0*p1];
% Set magnitude response
% at pole angle to be 1
z = exp(j*w0hat);
zvec = [z^2 z 1];
C = abs((denom*zvec') / (numer*zvec'));
freqz(C*numer, denom);

%%%% b.
%%%% Bandpass

%%%% c.
C4 = 261.626;
D4 = 293.665;
E4 = 329.628;
F4 = 349.228;
G4 = 391.995;
A4 = 440.000;
B4 = 493.883;
C5 = 523.251;
bpm = 60; % 300
beattime = 60/bpm;
fs = 8000;
 Ts = 1/fs;
N = beattime/Ts;
t = Ts : Ts : N*Ts;
f = [C4,D4,E4,F4,G4,A4,B4,C5];
vec = zeros(1, length(f)*N);
for i = 1:length(f)
    note = cos(2*pi*f(i)*t);
    vec((i-1)*N+1 : i*N) = note;
end
sound(vec, fs);


%%%% d.

N = beattime * fs;  % number of samples in one note
figure
spectrogram(vec, N, 100, N, fs, 'yaxis');
ylim( [0 0.7] );    % focus on 0 to 700 Hz


%%%% e.
y = filter(C*numer, denom, vec);

%%%% f.
pause(8);
sound(y,fs)
% it is more quiet and sounds cleaner

%%%%% g.
figure
spectrogram(y, N, 100, N, fs, 'yaxis');
ylim( [0 0.7] );    % focus on 0 to 700 Hz
% all the other frequencies that weren't the frequency of the note were
% filtered out leaving a much clearer power appearance of the specific
% frequency of the note
