%
% Calculate and Plot the Frequency Spectrum
% of a Signal with Optional Confidence Limits
%
%      [ h, f, sp ] = 
%        frspec( y, dt, pflag, clflag [, window ] )
%
% where   h      = plot handle
%         f      = frequency vector
%         sp     = frequency spectrum
%         y      = data time series
%         dt     = time interval of data
%         pflag  = plot if > 0
%         clflag = plot confidence limits if > 0
%         window = width of hanning window ( power
%                  of 2, default 8 )
%
% Bob Newell, February 1996
%
%----------------------------------------------------
%
function [ h, f, sp ] = ...
            frspec( y, dt, pflag, clflag, window )
%
% setup parameters
x = y(:);
n = length( x );
if nargin > 4, win = window; else win = 8; end; 
nfft = min( n, 2^win );
noverlap = 0;
fs = 2 * pi / dt;
% hanning window
window = 0.5 * ( 1 - cos( 2*pi*(1:nfft)'/(nfft+1) ) );
nwind = length( window );
k = fix( (n - noverlap)/(nwind-noverlap) );
index = 1:nwind;
kmu = k * norm( window )^2;
% for each window
pxx = zeros( nfft, 1 );
pxx2 = zeros( nfft, 1 );
a = [ (1:nfft)'/nfft ones(nfft,1) ];
for i = 1:k,
  xw = window .* ( x(index) - a * ( a \ x(index) ) );
  index = index + ( nwind - noverlap );
  xx = abs( fft( xw, nfft ) ) .^ 2;
  pxx = pxx + xx;
  pxx2 = pxx2 + abs( xx ) .^2;
end;
% select first half
if rem( nfft, 2 ),
  select = [ 1:(nfft+1)/2 ];
else,
  select = [ 1:nfft/2+1 ];
end;
pxx = pxx( select );
pxx2 = pxx2( select );
cpxx = zeros( size( pxx ) );
if k > 1,
  c = ( k .* pxx2 - abs( pxx ) .^ 2 ) ./ ( k - 1 );
  c = max( c, zeros( size( pxx ) ) );
  cpxx = sqrt( c );
end;
% confidence limits
ff = sqrt( 2 ) * erfinv( 0.95 );
pxx = pxx / kmu;
pxxc = ff .* cpxx / kmu;
sp = [ pxx pxxc ];
f = ( select - 1 )' * fs / nfft;
% plot results
if pflag,
  h = figure;
  if clflag,
    c = [ max( pxx - pxxc, 0 ) ( pxx + pxxc ) ];
    c = c .* ( c > 0 );
    plot( f, pxx, '-', f, c(:,1), '--', f, c(:,2), '--' )
  else,
    plot( f, pxx )
  end;
  set( gca, 'fontsize', 16 );
  xlabel( 'Frequency (radians/time)' );
else
  h = [];
end;
%
return;
%-----------------------------------------------------
% the end
