function show_lowpass
%
%   show_lowpass
%
%   Function to create the figures illustrating the low-pass filter and
%   convert to the frequency domain.

%Soupy Alexander, 3/31/04


figure
%Plot the filter weights
pl33=[  
  -0.00027 -0.00114 -0.00211 -0.00317 -0.00427 ...
  -0.00537 -0.00641 -0.00735 -0.00811 -0.00864 ...
  -0.00887 -0.00872 -0.00816 -0.00714 -0.00560 ...
  -0.00355 -0.00097  0.00213  0.00574  0.00980 ...
   0.01425  0.01902  0.02400  0.02911  0.03423 ...
   0.03923  0.04399  0.04842  0.05237  0.05576 ...
   0.05850  0.06051  0.06174  0.06215  ]';

% symmetric filter
pl33=[pl33(1:34); pl33(33:-1:1)];
plot([-33:33], pl33, '.')

%Set up the vector of zeros and a single one
theXs = [repmat(0, 1, 500) 1 repmat(0,1,500)];

%Run the data through the two filters
lp33Data = plfilt(theXs);

%Number of points for the FFT
nfft = 1024;

%If the data is longer than the filter, "wrap" it (repeat it) in sets equal
%to the number of FFT points.
if nfft < length(lp33Data);
    lp33Data = sum(buffer(lp33Data,nfft),2);
end
h33 = abs(fft(lp33Data,nfft));
w33 = 0:1/(nfft):1-1./nfft;

figure
loglog(w33, h33,'k')
l = line([1/24 1/24], [min(ylim) max(ylim)]);
set(l, 'color', 'k', 'linestyle', '-.');
t = text(1/24, 5, '24 Hours');
set(t, 'horiz', 'right', 'fontsize', 8)
l = line([1/12 1/12], [min(ylim) max(ylim)]);
set(l, 'color', 'k', 'linestyle', '-.');
t = text(1/12, 5, '12 Hours');
set(t, 'horiz', 'right', 'fontsize', 8)
l = line([1/6 1/6], [min(ylim) max(ylim)]);
set(l, 'color', 'k', 'linestyle', '-.');
t = text(1/6, 5, '6 Hours');
set(t, 'horiz', 'right', 'fontsize', 8)
xlabel('Frequency (cycles/h)')
ylabel(' ')

xlim([0.001 0.5])
ylim([0.0001 10])



theTime = [25.84 23.92 12.66 12.42 12.00];
theFreq = 1./theTime;
hitLocs = value2Index(w33, theFreq);
theDiffs = theFreq - w33(hitLocs);
theHValues = h33(hitLocs);