function  [sine1,sine2,sumsine,filtsum] = filtersinsefun(hh)

%%%  An impulse response is the input to this function.
%%%  
%%%  It plots the magnitude of the frequency response.
%%%
%%%  It samples and filters the sum of two sinewaves, and
%%%  outputs both sinewaves, their sum and the filtered sum.
%%%
%%%  It plots spectral lines corresponding to the two filtered sines.
%%%
%%%  It  plays in sequence: the first sinewave, the sum of the two
%%%  sinewawves, and the output of the filtered output.  
%%%
%%%  The frequencies f1 and f2 of the two sinewaves and
%%%  the sampling rate  fs  are specfied at the beginning of the script
%%%  

%%%   program parameters

fs = 2000 ;     % sampling frequency
T = 1.2 ;       % signal duration in seconds
f1 = 200 ;		% lower of the two frequencies
                % this is the frequency of the tone we wish to hear
f2 = 800 ;		% higher of the two frequencies, this is the frequency 
                % of the tone we wish to eliminate
Ts = 1/fs ;     % sampling interval
tt = 0:Ts:T ;   % time samples

dtf1 = 2*pi*f1*Ts ;
dtf2 = 2*pi*f2*Ts ;

disp('DT frequency of first sinusoid =')
disp(dtf1) ;
disp('DT frequency of second sinusoid =')
disp(dtf2) ;

sine1 = cos(2*pi*f1*tt) ;    % the first sinewave
sine2 = cos(2*pi*f2*tt) ;    % the second sinwave

sumsine = sine1 + sine2 ;    % the sum of sinewaves

filtsum = filter(hh,1,sumsine) ;  %the filtered sum

[HH,WW] = freqz(hh,1,200);    % the frequency response of the filter
hw = plot(WW,abs(HH),'k') ;
set(hw,'LineWidth',3)
hold on

n1 = (dtf1>WW) ;
n1 = sum(n1) ;
n2 = (dtf2>WW) ;
n2 = sum(n2) ;

hs = stem([dtf1 dtf2],[abs(HH(n1)),abs(HH(n2))]) ;  % plotting the spectral lines
set(hs,'LineWidth',3)
haz=gca ;
set(haz,'FontSize',18)
V=axis ;
axis([0 pi V(3) V(4)])
hx = xlabel('discrete-time frequency, 0 to pi') ;
set(hx,'FontSize',21)
hold off

pause

z=1;
switch z
case 1
	soundsc(sine1,fs)
	pause(1.5)   % needed for PC laptop, but not for my Mac
	soundsc(sumsine,fs)
	pause(1.5)
	soundsc(filtsum,fs)
case 2	%% in this case the output signals are upsampled before
	    %% playing.  This is intended to make permit a good reconstrution
		%% even when the sampling rate is little more than twice the highest freq.
	r = 2;
	xxx1 = resample(sine1,r,1);
	xxxx = resample(sumsine,r,1);
	yyyy = resample(filtsine,r,1);
	soundsc(xxx1,r*fs)
        pause(1.5)
	soundsc(xxxx,r*fs)
        pause(1.5)
	soundsc(yyyy,r*fs)
end




