% Computer homework #2
% T Kite March 1998

% Part 1 - subsampling

fs=8000; ts=1/fs;				% time vector
t=0:ts:ts*511;

a=sin(0.8*pi*t*fs)+sin(0.1*pi*t*fs);		% input signal
b=a(1:2:length(a));				% subsample
[ha,wa]=freqz(a.*hann(512),1,256,fs);		% compute frequency content
[hb,wb]=freqz(b.*hann(256),1,256,fs/2);

figure(1)
subplot(211); plot(wa,abs(ha))
xlabel 'Frequency (Hz)'; ylabel 'Magnitude'
text(1000,120,'Original frequency response')
subplot(212); plot(wb,abs(hb))
xlabel 'Frequency (Hz)'; ylabel 'Magnitude'
text(1000,60,'Subsampled frequency response')


% Part 2 - compensation for zero order hold

w=0:pi/128:255*pi/128;				% ideal response
fr=(w(1:129)+eps)./(2*sin((w(1:129)+eps)/2));	% eps prevents trouble at DC
fr(130:256)=conj(fr(128:-1:2));			% conjugate mirror
ir_ideal=fftshift(real(ifft(fr)));		% ideal impulse response

ir_trunc=ir_ideal(119:139);			% truncated impulse response
ir_trunc=ir_trunc/sum(ir_trunc);		% set DC gain to unity
ir_win=ir_trunc.*hamming(21)';			% ditto for windowed filter
ir_win=ir_win/sum(ir_win);

[h_trunc,w2]=freqz(ir_trunc,1,512,8000);	% truncated freq. response
h_trunc=abs(h_trunc);				% magnitude
h_win=abs(freqz(ir_win,1,512,8000));		% windowed freq. response
w3=linspace(0,pi,512)';
h_ideal=(w3+eps)./(2*sin((w3+eps)/2));		% ideal freq. resp. (512 pt)

figure(2)
subplot(211); plot(w2,h_trunc,w2,h_ideal,'--')
xlabel 'Frequency (Hz)'; ylabel 'Magnitude'
axis([0 4000 1 1.6]); legend('Truncated','Ideal',2)
subplot(212); plot(w2,h_win,w2,h_ideal,'--')
xlabel 'Frequency (Hz)'; ylabel 'Magnitude'
axis([0 4000 1 1.6]); legend('Windowed','Ideal',2)

h_trunc_dB=20*log10(h_trunc);			% freq. resp. in dB
h_win_dB=20*log10(h_win);
h_ideal_dB=20*log10(h_ideal);

figure(3)
subplot(211); act=plot(w2,h_trunc_dB-h_ideal_dB);
xlabel 'Frequency (Hz)'; ylabel 'Magnitude'
id=line([0 4000],[0 0]); set(id,'linestyle',':')
legend([act id],'Truncated','Ideal',3)
subplot(212); act=plot(w2,h_win_dB-h_ideal_dB);
xlabel 'Frequency (Hz)'; ylabel 'Magnitude'
id=line([0 4000],[0 0]); set(id,'linestyle',':')
legend([act id],'Windowed','Ideal',3)
