function bp = BPF(L, fc, fs, f)

% digital FIR bandpass filter system function, frequencies are in Hz.
% use: BPF(L, fc, fs, fstart, fstop, fpoints)
% L - number of samples in the impulse response sequence
% fc - center frequency
% fs - sample frequency
% fstart - fstop, range of frequencies for which to calculate HH
% fpoints - number of values in the HH array
%
% The -3dB bandwidth will be approximately 1/L (ie. fs/L Hz where fs is an actual
% sample frequency). Only positive frequencies are calculated and plotted.
%
w = 2*pi*f/fs; % scaled angular frequency array
wc = 2*pi*fc/fs; % scaled angular center frequency (for which abs(HH) has a maximum)
b = 2/L; % this adjusts the maximum value of the magnitude response to approximately 1

% The impulse response sequence is given as: h[n] = b*cos(wc*n).
% If the cosine function is written in terms of complex exponentials the
% frequency response function will consist of two parts, one for -wc and
% one for wc. Below, H1 is for wc and H2 for -wc. To obtain the complete
% response these two terms must be added.

H1 = 0.5*b*exp(j*0.5*(L - 1)*(wc - w)).*sin(0.5*L*(wc - w))./sin(0.5*(wc - w));

H2 = 0.5*b*exp(-j*0.5*(L - 1)*(wc + w)).*sin(0.5*L*(wc + w))./sin(0.5*(wc + w));

bp = H1 + H2;