function [y] = hfsgram20(x,Fs,preemp,dBmax,quefrencycutoff)
%	Calculate quality spectrogram for 16KHz data using
%       a 512-point FFT, and about a 0.03 sec cepstral
%       smoother, overlapping such that an analysis is done
%       every 10ms
%       Fix for top 32dB Resolution
%
%       Inputs
%       x - speech data 
%       Fs - sampling rate
%       preemp - preemphasis factor 0 = none, 1 = pure difference
%       dBmax  - how many top dB to display
%       quefrequencycutoff  - normalized quefrency for making wideband
%       spectrograms
%
%
%       H. Silverman 6-29-93
%       Copyright (c) 1993 -- LEMS, Brown University
%       Originally hacked by jtf for version 4.0
%       3/26/98 hacked by jea to implement cepstral smoothing 
%               and preemphasis
%       H. Silverman 12/15/20 -- make useful for 20kHz sampling

if nargin<2
  Fs = 20000.0;
  preemp = .5;
  dBmax=32;
  quefrencycutoff=.0015;
end

x = x(:); % make a column vector for ease later
nx=size(x,1)

%% do pre-emphasis now
x = x(2:nx) - preemp*x(1:nx-1);
nx=size(x,1);

nfft = 512;
nshift = round(.01*Fs); %% 10ms shift
noverlap = nfft-nshift;

%% Change the following line to change the window used.
win = hamming(nfft);

%% compute the segment offsets 
colindex = [1:nshift:nx-nfft];

%% Determine the Number of Spectra (Columns)
ncol = length(colindex);

%% preallocate the space for the dft's
y = zeros(nfft,ncol);

%% Major Loop on DFT

%% now stuff x into columns of y with the proper offset
%% should be able to do this with fancy indexing!

rowindex = [1:nfft]';
y(:) = x( rowindex(:,ones(1,ncol)) + colindex(ones(nfft,1),:)-1);

%% Apply the window to the array of offset signal segments.
y(:) = win(:,ones(1,ncol)).*y;

%% log mag take dft
y  = 20*log10(abs(fft(y,nfft)));

%% make a low-pass filter for cepstra
cepnfft = nfft/2
cepweighting = ones(cepnfft,1); %% cepstral weighting vector
quefrencies = [0:cepnfft-1]/Fs;  %% correpsonding quefrencies
cepweighting(find(quefrencies>quefrencycutoff),1)=0;
cepweighting(cepnfft:-1:cepnfft/2+2) = cepweighting(2:cepnfft/2); %% make symmetric

%% use window-design techniques to smooth this filter
cw = real(fft(cepweighting));
cw = cw([[cepnfft/2+1:cepnfft] [1:cepnfft/2]]) .* hamming(cepnfft);
cepweighting = abs(ifft(cw,nfft)); %% back to cepstral domain

yhat = real(ifft(y,nfft)); %% take to cepstra
yhat = cepweighting(:,ones(1,ncol)).*yhat; %% do weighting
y = real(fft(yhat)); %% back to log spectrum
y = y(1:nfft,:);

%% move down to max value then shift up dBmax dB
ymx = max(max(y));
y = y-ymx+dBmax;

fprintf(1,'max = %f min = %f\n',max(max(y)),min(min(y)));

%% zero out anything less than 0 dB
zerothese = find(y<0);
y(zerothese) = 0;

%% Set the Color Table

colormap(gray);

%%Image the spectrogram
imagesc(1000*colindex/Fs, [0 Fs/2], -y(1:nfft/2,:));
axis('xy');

%% Truncate axis

xlabel('time (ms)');
ylabel('freq (Hz)');