% 2007-10-09-freqresp.diary
%
% Passing a mixture of two sinusoid through a system with a known 
% frequency response (i.e. DTFT of its impulse response) modifies 
% the balance between the sinusoids in proportion to the frequency 
% response at those frequencies.

% 5 point moving average filter: frequency response
h = 1/5*ones(1,5);

% setup timebase
nn = -20:100;
% offset to get zero point
zz = 21;
% Mixture of two equal-amplitude sines at 0.1pi rad/samp and 0.5pi rad/samp
s1(zz+[0:100]) = sin(0.1*pi*[0:100]);
s2(zz+[0:100]) = sin(0.5*pi*[0:100]);
x = s1 + s2;

% Pass it through the system
y = conv(h,x);
% truncate y to be the same length as x
y = y(1:length(x));

% Plot results
subplot(411)
% The two component sinusoids
plot(nn,s1,nn,s2)
% Their sum (system input) and system output
subplot(412)
plot(nn,x,nn,y);
grid
axis([-20 100 -2 2])

% Compare to sinusoids scaled by freq resp
H = fft(h,256);
% with a 256pt fft, 2pi corresponds to 256, so 0.1pi corresponds to 256/20
H1 = H(1+round(0.1*(256/2)));
H2 = H(1+round(0.5*(256/2)));
% Individual sinusoids scaled and phase shifted by freq resp samples
s1b(zz+[0:100]) = abs(H1)*sin(0.1*pi*[0:100]+angle(H1));
s2b(zz+[0:100]) = abs(H2)*sin(0.5*pi*[0:100]+angle(H2));
% Plot the modified sinusoids
subplot(413)
plot(nn,s1b,nn,s2b)
% .. and compare their sum to the actual system output
subplot(414)
plot(nn,s1b+s2b,nn,y)
grid
axis([-20 100 -2 2])
% Pretty much fall fright on top of each other, except for a few
% points around zero (i.e. the transient response, which is limited to the
% impulse response length or 5 points).