function hw = HW(L, fc, fs)

% digital FIR bandpass filter system function
% use: HW(L (number of samples in the impulse response sequence), fc (resonance scaled frequency, f/fs),
% fs (sample frequency))
% 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.
%
f = linspace(0, 0.5, 500);
w = 2*pi*f;
wc = 2*pi*fc;
b = 1/fix((L + 1)/2); % this adjusts the maximum value of the magnitude response  to 1

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));

HH = H1 + H2;

subplot(2, 1, 1);

plot(f, abs(HH)), grid;

title('magnitude');

subplot(2, 1, 2);

plot(f, angle(HH)), grid;

title('phase');

f3dB = (find(abs(HH) >= 1/sqrt(2)) - 1)/499; % find frequency ranges where |HH| >= 0.7071.
% find(abs(HH) >= 1/sqrt(2)) only finds the indices pointing to the
% frequencies. The algorithm finds the scaled frequencies in the range 0 -
% 0.5. There are 500 sample points, ie there are then 499 spaces between
% the sample points.

f_Hz = fs*f3dB % print frquency values in Hz