function freqhist4(nf)

load /ligo/tmp2/bernard/matlab/results_20_points numFreq yy numBins nnr nni srate;

freq = (nf-1)*srate/2/(numFreq-1);
nframes = numBins*(numFreq-1)/25000;

%
% real part of fft
%q

N = sum(nnr(nf,2:end-1));
xx = yy(nf,:);
nn = nnr(nf,:);

% Calculate the variance and standard deviation sdr(nf) for each frequency nf.

[mur,sdr]=stat(xx,nn);

fprintf('mur = %f, sigmar = %f\n',  mur, sdr);

ny = gausscomp(mur, sdr, N, xx);
[mun,sdn]=stat(xx,ny);

subplot(2,1,1)

fprintf('mun = %f, sigman = %f\n',  mun, sdn);

doublehist(xx, nn, ny)

t = sprintf('Histograms at %f Hz (bin %i)', freq, nf);
title(t);
xlabel('Real part fft')
ylabel('Number of samples')
%
x1=(3*xx(1)-xx(2))/2;
x2=.9*max([nn,ny]);
x3=(8*xx(end)+xx(1))/9;
x4=.8*max([nn,ny]);
text(x1,x2,'red = data')
text(x1,x4,'blue = Gaussian fit')
text(x3,x2,['\mu = ',sprintf('%6.2f',mur)])
text(x3,x4,['\sigma = ',sprintf('%6.2f',sdr)])

%
% imaginary part of fft
%

N = sum(nni(nf,2:end-1));
nn = nni(nf,:);

% Calculate mean mui(nf) for each frequency nf.
% Calculate the variance and standard deviation sdr(nf) for each frequency nf.

[mui,sdi]=stat(xx,nn);

fprintf('mui = %f, sigmai = %f\n',  mui, sdi);

ny = gausscomp(mui, sdi, N, xx);
[mun,sdn]=stat(xx,ny);

fprintf('mun = %f, sigman = %f\n',  mun, sdn);


subplot(2,1,2)

doublehist(xx,nn,ny)
%hist(imag(bb(:,nf)),10)
xlabel('Imaginary part fft')
ylabel('Number of samples')
x1=(3*xx(1)-xx(2))/2;
x2=.9*max([nn,ny]);
x3=(8*xx(end)+xx(1))/9;
x4=.8*max([nn,ny]);
text(x1,x2,'red = data')
text(x1,x4,'blue = Gaussian fit')
text(x3,x2,['\mu = ',sprintf('%6.2f',mui)])
text(x3,x4,['\sigma = ',sprintf('%6.2f',sdi)])