%lowpassfilter.m
%Implements an ideal low pass digital filter as described in Section 5.6
% Impulse response h(n) is given in Section 5.6.3 and plotted in Fig. 5.24
% DTFT of the filter is given in Eq. 5.74
% Filter output y(t) is given in Eq. 5.58 (by convolution in time domain)

close all
b=pi/10; %bandwidth of the filter (for DTFT, pi is upper limit of freq.)
n=[-100:100]; %extent of impulse response to keep
h=(b/pi)*sinc(b.*n/pi); %impulse response
subplot(2,2,1); stem(n,h); %display impulse response
xlabel('discrete time, n')
title('Low-pass filter impulse response h[n]')

%----------------------------------------------------------------------
%--------------------Create input waveform x[n] (square wave)----------

%Let x(t) = square wave
% with fundamental frequency = 2/3 of the filter bandwidth.
%That way, the result of filtering should be a sinusoid at the
%fundamental frequency, since higher harmonics will be filtered out.

wx=(2/3)*b; %fundamental freq. of the input
perx=2*pi/wx; %period of input
perx2=perx/2; %half the period of the input
x1=ones(1,perx2); %square wave = 1 for half the period
x0=-1*ones(1,perx2);%square wave = -1 for half the period
x=[x1 x0]; %one period of the input
ln=length(n); %number of time samples
while (length(x) < ln); %generate about as many samples of the input
%as we already have of the filter impulse response
x=[x x]; %replicate periods to build the input square wave
end
lx=length(x); lx2=lx/2;ax=[-lx2:lx2-1]; %center the time axis
subplot(2,2,2); stem(ax,x); %plot the input square wave
xlabel('discrete time, n')
title('input square wave')

%-----------------------------------------------------------------------
%-------Now compute the output using convolution in the time domain-----

y=conv(x,h); %compute the output
ly=length(y); ly2=ly/2;ay=[-ly2:ly2-1]; %center the time axis
subplot(2,2,3); stem(ay,y); %stem plot the output
xlabel('discrete time, n')
title('Stem plot of filtered output y[n]')

subplot(2,2,4);plot(ay,y); %smoothly plot the output
xlabel('discrete time, n')
title('Smooth plot of filtered output y[n]')