%
% compare(f, flag)
%
% Draws a pair of histograms showing differences before and after line
% removal at a particular frequency
% 
% If flag = 0 or is left undefined, f is interpreted as a frequency
% number in the range 0 to numFreq
%
% If flag = 1, f is interpreted as frequency between 0 and the Nyquist
% frequency.
%
function comp40LR(f, flag)
%load /ligo/tmp2/bernard/matlab/19nov94LR2/results_before
load /ligo/tmp2/bernard/matlab/results_1000 yy nnr nni srate fftPad numFreq;

if (nargin < 2)
    flag = 0;
end;

if (flag == 0)
    if (f < 1 | f > numFreq)
       error(sprintf('f must lie between 1 and %d', numFreq));
    end;
    nf = f;
else
    if (f < 0 | f > srate/2)
       error(sprintf('f must lie between 0 and %f', srate/2));
    end;
    nf = round(fftPad*f/srate) + 1;
end;

freq = (nf-1)*srate/fftPad;

x = yy(nf,:);
nr1 = nnr(nf,:);
ni1 = nni(nf,:);

%load /ligo/tmp2/bernard/matlab/19nov94LR2/results_after
load /ligo/tmp2/bernard/matlab/results_afterLR yy nnr nni srate fftPad blockTotal numPoints numFreq;

nr2 = nnr(nf,:);
ni2 = nni(nf,:);

subplot(2,1,1);
doublehist(x, nr1, nr2);

title(sprintf('Histograms at %.2f Hz (f_n=%d) from %d points', freq, nf, blockTotal*numPoints));

height = max([nnr(nf,:)]);
axis([1.05*yy(nf,1) 1.05*yy(nf,end) 0 height]);
x1 = yy(nf,1) + 0.05*(yy(nf,end)-yy(nf,1));
x2 = 0.9*height;
x4 = 0.8*height;

text(x1, x2, 'red = before LR');
text(x1, x4, 'blue = after LR');

subplot(2,1,2);

doublehist(x, ni1, ni2);

height = max([nni(nf,:)]);
axis([1.05*yy(nf,1) 1.05*yy(nf,end) 0 height]);
x1 = yy(nf,1) + 0.05*(yy(nf,end)-yy(nf,1));
x2 = 0.9*height;
x4 = 0.8*height;

text(x1, x2, 'red = before LR');
text(x1, x4, 'blue = after LR');