%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "Insert title here" 
% Written for GNU Octave
% 
% INSTITUTE OF TELECOMMUNICATIONS  
% University of Stuttgart 
% www.inue-uni-stuttgart.de 
% author: Your Name   
% date: DD.MM.YYYY   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This example script just prints the values of the given parameters.
% Make sure that the function prototype matches the prototype defined
% in the XML config file. (Parameter names are arbitrary but should be the same
% for better readability)

function [] = demo_2_fun(select1, enum1, window_size, filename)

% define formats
landscape = "-S930,350";
portrait  = "-S640,480";

% select format
output_format = landscape;

% Create an invisible figure.
fh = figure(1);
set(fh, "visible", "off");

% define style
set(0, "defaultlinelinewidth", 2);
set(0, "defaultaxesfontsize", 8);
set(0, "defaultaxesfontname", "Helvetica");
set(0, "defaulttextfontsize", 8);
set(0, "defaulttextfontname", "Helvetica");

% Plot into the figure


if (strcmpi(enum1,'psd'))
  psd = 1;
  spectrogram = 0;
end
if (strcmpi(enum1,'spectrogram'))
  psd = 0;
  spectrogram = 1;
end
if (strcmpi(enum1,'both'))
  psd = 1;
  spectrogram = 1;
end


if (select1 == 1)
  wave_file='wlan.Wfm.bin';
elseif (select1 == 2)
  wave_file='bluetooth.Wfm.bin';
elseif (select1 == 3)
  wave_file='micro.Wfm.bin';
elseif (select1 == 4)
  wave_file='bluetooth_wlan.Wfm.bin';
elseif (select1 == 5)
  wave_file='microwave_wlan.Wfm.bin';
elseif (select1 == 6)
  wave_file='all.Wfm.bin';
elseif (select1 == 7)
  wave_file='wlan2.Wfm.bin';
elseif (select1 == 8)
  wave_file='bluetooth_wlan2.Wfm.bin';  
end

v=load_bin(wave_file);
if (select1 == 1)
  v=v.*4;
end



if (psd == 1 && spectrogram == 1)
  sub_ax1 = axes('Position', [0.07 0.16 0.39 0.72])
  %sub1=subplot(1,2,1)
  %set (sub1, 'Position', [0 2 0.7 0.7] );
end

if (psd == 1)
  pwelch(v,1000);
  XTick = linspace(0,0.5,11);
  XTickLabel = linspace(2.4,2.5,11);
  ylabel('measured spectral power in dB')
  xlabel('frequency in GHz');
%  if (select1 == 1)
%    ylim([0.00000001 1]);
%    YTick = logspace(-8,0,8);
%    YTickLabel = -80:10:0;
%  else
    ylim([0.0000001 1]);
    YTick = logspace(-7,0,8);
    YTickLabel = -70:10:0;
%  end
  title('power spectral density');
  set(gca,'yscale','log','XTick',XTick,'XTickLabel',XTickLabel,'YTick',YTick,'YTickLabel',YTickLabel); 
end  

if (psd == 1 && spectrogram == 1)
  sub_ax2 = axes('Position', [0.56 0.16 0.43 0.72])
  %subplot(1,2,2)
elseif (psd == 0 && spectrogram == 1)
  sub_ax2 = axes('Position', [0.1 0.16 0.83 0.72])
end
if (spectrogram == 1)
  specgram(v,window_size,2,hanning(window_size));
  title('spectrogram');
  ylabel('frequency in GHz');
  xlabel('time in ms');
  XTick = linspace(0,1e6,11);
  XTickLabel = linspace(0,10,11);
  YTick = linspace(0,0.99,11);
  YTickLabel = linspace(2.4,2.5,11);
  set(gca,'XTick',XTick,'XTickLabel',XTickLabel,'YTick',YTick,'YTickLabel',YTickLabel);
  
  if (psd == 0 && spectrogram == 1)
    cb = colorbar;
    zlab = get(cb,'ylabel');
    set(zlab,'String',['power' blanks(1) 'density' blanks(1) 'in' blanks(1) 'dB']); 
  end
end


% Create an image file. This png image is displayed on the webpage. Its size 
% must be 930x350.
print(filename, "-dpng", output_format);

end









%% needed functions..


function ret = roundn (x, n = 0)

  if (mod (x, 1) != 0)

    ret = round (10^abs (n) * x) / (10^abs (n));
    
  else
  
    ret = round (x / 10^abs (n)) * 10 ^ abs (n);
  
  end  

end




function [v]=load_bin(name)
    fid=fopen(name,'rb');
    fread(fid, 2,'single=>double');
    v=fread(fid, inf,'single=>double');
    fclose(fid);
end



function [S_r, f_r, t_r] = specgram(x, n = min(256, length(x)), Fs = 2, window = hanning(n), overlap = ceil(length(window)/2))

  if nargin < 1 || nargin > 5
    print_usage;
  ## make sure x is a vector
  elseif columns(x) != 1 && rows(x) != 1
    error ("specgram data must be a vector");
  endif
  if columns(x) != 1, x = x'; endif

  ## if only the window length is given, generate hanning window
  if length(window) == 1, window = hanning(window); endif

  ## should be extended to accept a vector of frequencies at which to
  ## evaluate the Fourier transform (via filterbank or chirp
  ## z-transform)
  if length(n)>1,
    error("specgram doesn't handle frequency vectors yet");
  endif

  if (length (x) <= length (window))
    error ("specgram: segment length must be less than the size of X");
  endif

  ## compute window offsets
  win_size = length(window);
  if (win_size > n)
    n = win_size;
    warning ("specgram fft size adjusted to %d", n);
  endif
  step = win_size - overlap;

  ## build matrix of windowed data slices
  offset = [ 1 : step : length(x)-win_size ];
  S = zeros (n, length(offset));
  for i=1:length(offset)
    S(1:win_size, i) = x(offset(i):offset(i)+win_size-1) .* window;
  endfor

  ## compute Fourier transform
  S = fft (S);

  ## extract the positive frequency components
  if rem(n,2)==1
    ret_n = (n+1)/2;
  else
    ret_n = n/2;
  endif
  S = S(1:ret_n, :);

  f = [0:ret_n-1]*Fs/n;
  t = offset/Fs;
  if nargout==0
    imagesc(t, f, 20*log10(abs(S)));
    set (gca (), "ydir", "normal");
    xlabel ("Time")
    ylabel ("Frequency")
  endif
  if nargout>0, S_r = S; endif
  if nargout>1, f_r = f; endif
  if nargout>2, t_r = t; endif

end

function varargout = pwelch(x,varargin)

  %%
  %% COMPATIBILITY LEVEL
  %% Argument positions and defaults depend on compatibility level selected
  %% by calling pwelch without arguments or with a single string argument.
  %%   native:      compatib=1; prev_compat=pwelch(); prev_compat=pwelch([]);
  %%   matlab R11:  compatib=2; prev_compat=pwelch('R11-');
  %%   matlab R12:  compatib=3; prev_compat=pwelch('R12+');
  %%   spectrum.welch defaults:  compatib=4; prev_compat=pwelch('psd');
  %% In each case, the returned value is the PREVIOUS compatibility string.
  %%
  compat_str = {[]; 'R11-'; 'R12+'; 'psd'};
  persistent compatib;
  if ( isempty(compatib) || compatib<=0 || compatib>4 )
    %% legal values are 1, 2, 3, 4

    compatib = 1;
  endif
  if ( nargin <= 0 )
    error( 'pwelch: Need at least 1 arg. Use "help pwelch".' );
  elseif ( nargin==1 && (ischar(x) || isempty(x)) )
    varargout{1} = compat_str{compatib};
    if ( isempty(x) ) # native
      compatib = 1;
    elseif ( strcmp(x,'R11-') )
      compatib = 2;
    elseif ( strcmp(x,'R12+') )
      compatib = 3;
    elseif ( strcmp(x,'psd') )
      compatib = 4;
    else
      error( 'pwelch: compatibility arg must be empty, R11-, R12+ or psd' );
    endif
    %% return
  %%
  %% Check fixed argument
  elseif ( isempty(x) || ~isvector(x) )
    error( 'pwelch: arg 1 (x) must be vector.' );
  else
    %%  force x to be COLUMN vector
    if ( size(x,1)==1 )
      x=x(:);
    endif
    %%
    %% Look through all args to check if  cross PSD, transfer function or
    %% coherence is required.  If yes, the second arg is data vector "y".
    arg2_is_y = 0;
    x_len = length(x);
    nvarargin = length(varargin);
    for iarg=1:nvarargin
      arg = varargin{iarg};
      if ( ~isempty(arg) && ischar(arg) && ...
           ( strcmp(arg,'cross') || strcmp(arg,'trans') || ...
             strcmp(arg,'coher') || strcmp(arg,'ypower') ))
        %% OK. Need "y". Grab it from 2nd arg.
        arg = varargin{1};
        if ( nargin<2 || isempty(arg) || ~isvector(arg) || length(arg)~=x_len )
          error( 'pwelch: arg 2 (y) must be vector, same length as x.' );
        endif
        %% force  COLUMN vector
        y = varargin{1}(:);
        arg2_is_y = 1;
        break;
      endif
    endfor
    %%
    %% COMPATIBILITY
    %% To select default argument values, "compatib" is used as an array index.
    %% Index values are   1=native,  2=R11,  3=R12,  4=spectrum.welch
    %%
    %%  argument positions:
    %%  arg_posn = varargin index of window, overlap, Nfft, Fs and conf
    %%             args respectively, a value of zero ==>> arg does not exist
    arg_posn = [1 2 3 4 5;  # native
                3 4 1 2 5;  # Matlab R11- pwelch
                1 2 3 4 0;  # Matlab R12+ pwelch
                1 2 3 4 5]; # spectrum.welch defaults
    arg_posn  = arg_posn(compatib,:) + arg2_is_y;
    %%
    %%  SPECIFY SOME DEFAULT VALUES for (not all) optional arguments
    %%  Use compatib as array index.
    %%  Fs = sampling frequency
    Fs        = [ 1.0 2*pi 2*pi 2*pi ];
    Fs        = Fs(compatib);
    %%  plot_type: 1='plot'|'squared'; 5='db'|'dB'
    plot_type = [ 1 5 5 5 ];
    plot_type = plot_type(compatib);
    %%  rm_mean: 3='long-mean'; 0='no-strip'|'none'
    rm_mean   = [ 3 0 0 0 ];
    rm_mean   = rm_mean(compatib);
    %% use max_overlap=x_len-1 because seg_len is not available yet
    %% units of overlap are different for each version:
    %%    fraction, samples, or percent
    max_overlap = [ 0.95 x_len-1 x_len-1 95];
    max_overlap = max_overlap(compatib);
    %% default confidence interval
    %%  if there are more than 2 return values and if there is a "conf" arg
    conf      = 0.95 * (nargout>2) * (arg_posn(5)>0);
    %%
    is_win    = 0;    # =0 means valid window arg is not provided yet
    Nfft      = [];   # default depends on segment length
    overlap   = [];   # WARNING: units can be #samples, fraction or percentage
    range     = ~isreal(x) || ( arg2_is_y && ~isreal(y) );
    is_sloppy = 0;
    n_results = 0;
    do_power  = 0;
    do_cross  = 0;
    do_trans  = 0;
    do_coher  = 0;
    do_ypower = 0;
  %%
  %%  DECODE AND CHECK OPTIONAL ARGUMENTS
    end_numeric_args = 0;
    for iarg = 1+arg2_is_y:nvarargin
      arg = varargin{iarg};
      if ( ischar(arg) )
        %% first string arg ==> no more numeric args
        %% non-string args cannot follow a string arg
        end_numeric_args = 1;
        %%
        %% decode control-string arguments
        if ( strcmp(arg,'sloppy') )
          is_sloppy = ~is_win || is_win==1;
        elseif ( strcmp(arg,'plot') || strcmp(arg,'squared') )
          plot_type = 1;
        elseif ( strcmp(arg,'semilogx') )
          plot_type = 2;
        elseif ( strcmp(arg,'semilogy') )
          plot_type = 3;
        elseif ( strcmp(arg,'loglog') )
          plot_type = 4;
        elseif ( strcmp(arg,'db') || strcmp(arg,'dB') )
          plot_type = 5;
        elseif ( strcmp(arg,'half') || strcmp(arg,'onesided') )
          range = 0;
        elseif ( strcmp(arg,'whole') || strcmp(arg,'twosided') )
          range = 1;
        elseif ( strcmp(arg,'shift') || strcmp(arg,'centerdc') )
          range = 2;
        elseif ( strcmp(arg,'long-mean') )
          rm_mean = 3;
        elseif ( strcmp(arg,'linear') )
          rm_mean = 2;
        elseif ( strcmp(arg,'short') || strcmp(arg,'mean') )
          rm_mean = 1;
        elseif ( strcmp(arg,'no-strip') || strcmp(arg,'none') )
          rm_mean = 0;
        elseif ( strcmp(arg, 'power' ) )
          if ( ~do_power )
            n_results = n_results+1;
            do_power = n_results;
          endif
        elseif ( strcmp(arg, 'cross' ) )
          if ( ~do_cross )
            n_results = n_results+1;
            do_cross = n_results;
          endif
        elseif ( strcmp(arg, 'trans' ) )
          if ( ~do_trans )
            n_results = n_results+1;
            do_trans = n_results;
          endif
        elseif ( strcmp(arg, 'coher' ) )
          if ( ~do_coher )
            n_results = n_results+1;
            do_coher = n_results;
          endif
        elseif ( strcmp(arg, 'ypower' ) )
          if ( ~do_ypower )
            n_results = n_results+1;
            do_ypower = n_results;
          endif
        else
          error( 'pwelch: string arg %d illegal value: %s', iarg+1, arg );
        endif
        %% end of processing string args
        %%
      elseif ( end_numeric_args )
        if ( ~isempty(arg) )
          %% found non-string arg after a string arg ... oops
          error( 'pwelch: control arg must be string' );
        endif
      %%
      %% first 4 optional arguments are numeric -- in fixed order
      %%
      %% deal with "Fs" and "conf" first because empty arg is a special default
      %% -- "Fs" arg -- sampling frequency
      elseif ( iarg == arg_posn(4) )
        if ( isempty(arg) )
          Fs = 1;
        elseif ( ~isscalar(arg) || ~isreal(arg) || arg<0 )
          error( 'pwelch: arg %d (Fs) must be real scalar >0', iarg+1 );
        else
          Fs = arg;
        endif
      %%
      %%  -- "conf" arg -- confidence level
      %%    guard against the "it cannot happen" iarg==0
      elseif ( arg_posn(5) && iarg == arg_posn(5) )
        if ( isempty(arg) )
          conf = 0.95;
        elseif ( ~isscalar(arg) || ~isreal(arg) || arg < 0.0 || arg >= 1.0 )
          error( 'pwelch: arg %d (conf) must be real scalar, >=0, <1',iarg+1 );
        else
          conf = arg;
        endif
      %%
      %% skip all empty args from this point onward
      elseif ( isempty(arg) )
        1;
      %%
      %%  -- "window" arg -- window function
      elseif ( iarg == arg_posn(1) )
        if ( isscalar(arg) )
          is_win = 1;
        elseif ( isvector(arg) )
          is_win = length(arg);
          if ( size(arg,2)>1 )  # vector must be COLUMN vector
            arg = arg(:);
          endif
        else
          is_win = 0;
        endif
        if ( ~is_win )
          error( 'pwelch: arg %d (window) must be scalar or vector', iarg+1 );
        elseif ( is_win==1 && ( ~isreal(arg) || fix(arg)~=arg || arg<=3 ) )
          error( 'pwelch: arg %d (window) must be integer >3', iarg+1 );
        elseif ( is_win>1 && ( ~isreal(arg) || any(arg<0) ) )
          error( 'pwelch: arg %d (window) vector must be real and >=0',iarg+1);
        endif
        window = arg;
        is_sloppy = 0;
      %%
      %% -- "overlap" arg -- segment overlap
      elseif ( iarg == arg_posn(2) )
        if (~isscalar(arg) || ~isreal(arg) || arg<0 || arg>max_overlap )
          error( 'pwelch: arg %d (overlap) must be real from 0 to %f', ...
                 iarg+1, max_overlap );
        endif
        overlap = arg;
      %%
      %% -- "Nfft" arg -- FFT length
      elseif ( iarg == arg_posn(3) )
        if ( ~isscalar(arg) || ~isreal(arg) || fix(arg)~=arg || arg<0 )
          error( 'pwelch: arg %d (Nfft) must be integer >=0', iarg+1 );
        endif
        Nfft = arg;
      %%
      else
        error( 'pwelch: arg %d  must be string', iarg+1 );
      endif
    endfor
    if ( conf>0 && (n_results && ~do_power ) )
      error('pwelch: can give confidence interval for x power spectrum only' );
    endif
    %%
    %% end DECODE AND CHECK OPTIONAL ARGUMENTS.
    %%
    %% SETUP REMAINING PARAMETERS
    %% default action is to calculate power spectrum only
    if ( ~n_results )
      n_results = 1;
      do_power = 1;
    endif
    need_Pxx = do_power || do_trans || do_coher;
    need_Pxy = do_cross || do_trans || do_coher;
    need_Pyy = do_coher || do_ypower;
    log_two = log(2);
    nearly_one = 0.99999999999;
    %%
    %% compatibility-options
    %% provides exact compatibility with Matlab R11 or R12
    %%
    %% Matlab R11 compatibility
    if ( compatib==2 )
      if ( isempty(Nfft) )
        Nfft = min( 256, x_len );
      endif
      if ( is_win > 1 )
        seg_len = min( length(window), Nfft );
        window = window(1:seg_len);
      else
        if ( is_win )
          %% window arg is scalar
          seg_len = window;
        else
          seg_len = Nfft;
        endif
        %% make Hann window (don't depend on sigproc)
        xx = seg_len - 1;
        window = 0.5 - 0.5 * cos( (2*pi/xx)*[0:xx].' );
      endif
    %%
    %% Matlab R12 compatibility
    elseif ( compatib==3 )
      if ( is_win > 1 )
        %% window arg provides window function
        seg_len = length(window);
      else
        %% window arg does not provide window function; use Hamming
        if ( is_win )
          %% window arg is scalar
          seg_len = window;
        else
          %% window arg not available; use R12 default, 8 windows
          %% ignore overlap arg; use overlap=50% -- only choice that makes sense
          %% this is the magic formula for 8 segments with 50% overlap
          seg_len = fix( (x_len-3)*2/9 );
        endif
        %% make Hamming window (don't depend on sigproc)
        xx = seg_len - 1;
        window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' );
      endif
      if ( isempty(Nfft) )
        Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) );
      endif
    %%
    %% Matlab R14 psd(spectrum.welch) defaults
    elseif ( compatib==4 )
      if ( is_win > 1 )
        %% window arg provides window function
        seg_len = length(window);
      else
        %% window arg does not provide window function; use Hamming
        if ( is_win )
          %% window arg is scalar
          seg_len = window;
        else
          %% window arg not available; use default seg_len = 64
          seg_len = 64;
        endif
        %% make Hamming window (don't depend on sigproc)
        xx = seg_len - 1;
        window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' );
      endif
      %% Now we know segment length,
      %% so we can set default overlap as number of samples
      if ( ~isempty(overlap) )
        overlap = fix(seg_len * overlap / 100 );
      endif
      if ( isempty(Nfft) )
        Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) );
      endif
    %%
    %% default compatibility level
    else # if ( compatib==1 )
      %% calculate/adjust segment length, window function
      if ( is_win > 1 )
        %% window arg provides window function
        seg_len = length(window);
      else
        %% window arg does not provide window function; use Hamming
        if ( is_win )       # window arg is scalar
          seg_len = window;
        else
          %% window arg not available; use default length:
          %% = sqrt(length(x)) rounded up to nearest integer power of 2
          if ( isempty(overlap) )
            overlap=0.5;
          endif
          seg_len = 2 ^ ceil( log(sqrt(x_len/(1-overlap)))*nearly_one/log_two );
        endif
        %% make Hamming window (don't depend on sigproc)
        xx = seg_len - 1;
        window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' );
      endif
      %% Now we know segment length,
      %% so we can set default overlap as number of samples
      if ( ~isempty(overlap) )
        overlap = fix(seg_len * overlap);
      endif
      %%
      %% calculate FFT length
      if ( isempty(Nfft) )
        Nfft = seg_len;
      endif
      if ( is_sloppy )
        Nfft = 2 ^ ceil( log(Nfft) * nearly_one / log_two );
      endif
    endif
    %% end of compatibility options
    %%
    %% minimum FFT length is seg_len
    Nfft = max( Nfft, seg_len );
    %% Mean square of window is required for normalizing PSD amplitude.
    win_meansq = (window.' * window) / seg_len;
    %%
    %% Set default or check overlap.
    if ( isempty(overlap) )
      overlap = fix(seg_len /2);
    elseif ( overlap >= seg_len )
      error( 'pwelch: arg (overlap=%d) too big. Must be <length(window)=%d',...
             overlap, seg_len );
    endif
    %%
    %% Pad data with zeros if shorter than segment. This should not happen.
    if ( x_len < seg_len )
      x = [x; zeros(seg_len-x_len,1)];
      if ( arg2_is_y )
        y = [y; zeros(seg_len-x_len,1)];
      endif
      x_len = seg_len;
    endif
    %% end SETUP REMAINING PARAMETERS
    %%
    %%
    %% MAIN CALCULATIONS
    %% Remove mean from the data
    if ( rm_mean == 3 )
      n_ffts = max( 0, fix( (x_len-seg_len)/(seg_len-overlap) ) ) + 1;
      x_len  = min( x_len, (seg_len-overlap)*(n_ffts-1)+seg_len );
      if ( need_Pxx || need_Pxy )
        x = x - sum( x(1:x_len) ) / x_len;
      endif
      if ( arg2_is_y || need_Pxy)
        y = y - sum( y(1:x_len) ) / x_len;
      endif
    endif
    %%
    %% Calculate and accumulate periodograms
    %%   xx and yy are padded data segments
    %%   Pxx, Pyy, Pyy are periodogram sums, Vxx is for confidence interval
    xx = zeros(Nfft,1);
    yy = xx;
    Pxx = xx;
    Pxy = xx;
    Pyy = xx;
    if ( conf>0 )
      Vxx = xx;
    else
      Vxx = [];
    endif
    n_ffts = 0;
    for start_seg = [1:seg_len-overlap:x_len-seg_len+1]
      end_seg = start_seg+seg_len-1;
      %% Don't truncate/remove the zero padding in xx and yy
      if ( need_Pxx || need_Pxy )
        if ( rm_mean==1 ) # remove mean from segment
          xx(1:seg_len) = window .* ( ...
            x(start_seg:end_seg) - sum(x(start_seg:end_seg)) / seg_len);
        elseif ( rm_mean == 2 ) # remove linear trend from segment
          xx(1:seg_len) = window .* detrend( x(start_seg:end_seg) );
        else # rm_mean==0 or 3
          xx(1:seg_len) = window .* x(start_seg:end_seg);
        endif
        fft_x = fft(xx);
      endif
      if ( need_Pxy || need_Pyy )
        if ( rm_mean==1 ) # remove mean from segment
          yy(1:seg_len) = window .* ( ...
            y(start_seg:end_seg) - sum(y(start_seg:end_seg)) / seg_len);
        elseif ( rm_mean == 2 ) # remove linear trend from segment
          yy(1:seg_len) = window .* detrend( y(start_seg:end_seg) );
        else # rm_mean==0 or 3
          yy(1:seg_len) = window .* y(start_seg:end_seg);
        endif
        fft_y = fft(yy);
      endif
      if ( need_Pxx )
        %% force Pxx to be real; pgram = periodogram
        pgram = real(fft_x .* conj(fft_x));
        Pxx = Pxx + pgram;
        %% sum of squared periodograms is required for confidence interval
        if ( conf>0 )
          Vxx = Vxx + pgram .^2;
        endif
      endif
      if ( need_Pxy )
        %% Pxy (cross power spectrum) is complex. Do not force to be real.
        Pxy = Pxy + fft_y .* conj(fft_x);
      endif
      if ( need_Pyy )
        %% force Pyy to be real
        Pyy = Pyy + real(fft_y .* conj(fft_y));
      endif
      n_ffts = n_ffts +1;
    endfor
    %%
    %% Calculate confidence interval
    %%    -- incorrectly assumes that the periodogram has Gaussian probability
    %%       distribution (actually, it has a single-sided (e.g. exponential)
    %%       distribution.
    %% Sample variance of periodograms is (Vxx-Pxx.^2/n_ffts)/(n_ffts-1).
    %%    This method of calculating variance is more susceptible to round-off
    %%  error, but is quicker, and for double-precision arithmetic and the
    %%  inherently noisy periodogram (variance==mean^2), it should be OK.
    if ( conf>0 && need_Pxx )
      if ( n_ffts<2 )
        Vxx = zeros(Nfft,1);
      else
        %% Should use student distribution here (for unknown variance), but tinv
        %% is not a core Matlab function (is in statistics toolbox. Grrr)
        Vxx = (erfinv(conf)*sqrt(2*n_ffts/(n_ffts-1))) * sqrt(Vxx-Pxx.^2/n_ffts);
      endif
    endif
    %%
    %% Convert two-sided spectra to one-sided spectra (if range == 0).
    %% For one-sided spectra, contributions from negative frequencies are added
    %% to the positive side of the spectrum -- but not at zero or Nyquist
    %% (half sampling) frequencies.  This keeps power equal in time and spectral
    %% domains, as required by Parseval theorem.
    %%
    if ( range == 0 )
      if ( ~ rem(Nfft,2) )    # one-sided, Nfft is even
        psd_len = Nfft/2+1;
        if ( need_Pxx )
          Pxx = Pxx(1:psd_len) + [0; Pxx(Nfft:-1:psd_len+1); 0];
          if ( conf>0 )
            Vxx = Vxx(1:psd_len) + [0; Vxx(Nfft:-1:psd_len+1); 0];
          endif
        endif
        if ( need_Pxy )
          Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1); 0]);
        endif
        if ( need_Pyy )
          Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1); 0];
        endif
      else                    # one-sided, Nfft is odd
        psd_len = (Nfft+1)/2;
        if ( need_Pxx )
          Pxx = Pxx(1:psd_len) + [0; Pxx(Nfft:-1:psd_len+1)];
          if ( conf>0 )
            Vxx = Vxx(1:psd_len) + [0; Vxx(Nfft:-1:psd_len+1)];
          endif
        endif
        if ( need_Pxy )
          Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1)]);
        endif
        if ( need_Pyy )
          Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1)];
        endif
      endif
    else                      # two-sided (and shifted)
      psd_len = Nfft;
    endif
    %% end MAIN CALCULATIONS
    %%
    %% SCALING AND OUTPUT
    %% Put all results in matrix, one row per spectrum
    %%   Pxx, Pxy, Pyy are sums of periodograms, so "n_ffts"
    %%   in the scale factor converts them into averages
    spectra    = zeros(psd_len,n_results);
    spect_type = zeros(n_results,1);
    scale = n_ffts * seg_len * Fs * win_meansq;
    if ( do_power )
      spectra(:,do_power) = Pxx / scale;
      spect_type(do_power) = 1;
      if ( conf>0 )
        Vxx = [Pxx-Vxx Pxx+Vxx]/scale;
      endif
    endif
    if ( do_cross )
      spectra(:,do_cross) = Pxy / scale;
      spect_type(do_cross) = 2;
    endif
    if ( do_trans )
      spectra(:,do_trans) = Pxy ./ Pxx;
      spect_type(do_trans) = 3;
    endif
    if ( do_coher )
      %% force coherence to be real
      spectra(:,do_coher) = real(Pxy .* conj(Pxy)) ./ Pxx ./ Pyy;
      spect_type(do_coher) = 4;
    endif
    if ( do_ypower )
      spectra(:,do_ypower) = Pyy / scale;
      spect_type(do_ypower) = 5;
    endif
    freq = [0:psd_len-1].' * ( Fs / Nfft );
    %%
    %% range='shift': Shift zero-frequency to the middle
    if ( range == 2 )
      len2 = fix((Nfft+1)/2);
      spectra = [ spectra(len2+1:Nfft,:); spectra(1:len2,:)];
      freq    = [ freq(len2+1:Nfft)-Fs; freq(1:len2)];
      if ( conf>0 )
        Vxx = [ Vxx(len2+1:Nfft,:); Vxx(1:len2,:)];
      endif
    endif
    %%
    %%  RETURN RESULTS or PLOT
    if ( nargout>=2 && conf>0 )
      varargout{2} = Vxx;
    endif
    if ( nargout>=(2+(conf>0)) )
      %% frequency is 2nd or 3rd return value,
      %% depends on if 2nd is confidence interval
      varargout{2+(conf>0)} = freq;
    endif
    if ( nargout>=1 )
      varargout{1} = spectra;
    else
      %%
      %% Plot the spectra if there are no return variables.
      plot_title=['power spectrum x ';
                  'cross spectrum   ';
                  'transfer function';
                  'coherence        ';
                  'power spectrum y ' ];
      for ii = 1: n_results
        if ( conf>0 && spect_type(ii)==1 )
          Vxxxx = Vxx;
        else
          Vxxxx = [];
        endif
        if ( n_results > 1 )
          figure();
        endif
        if ( plot_type == 1 )
          plot(freq,[abs(spectra(:,ii)) Vxxxx]);
        elseif ( plot_type == 2 )
          semilogx(freq,[abs(spectra(:,ii)) Vxxxx]);
        elseif ( plot_type == 3 )
          semilogy(freq,[abs(spectra(:,ii)) Vxxxx]);
        elseif ( plot_type == 4 )
          loglog(freq,[abs(spectra(:,ii)) Vxxxx]);
        elseif ( plot_type == 5 )  # db
          ylabel( 'amplitude (dB)' );
          plot(freq,[10*log10(abs(spectra(:,ii))) 10*log10(abs(Vxxxx))]);
        endif
        title( char(plot_title(spect_type(ii),:)) );
        ylabel( 'amplitude' );
        %% Plot phase of cross spectrum and transfer function
        if ( spect_type(ii)==2 || spect_type(ii)==3 )
          figure();
          if ( plot_type==2 || plot_type==4 )
            semilogx(freq,180/pi*angle(spectra(:,ii)));
          else
            plot(freq,180/pi*angle(spectra(:,ii)));
          endif
          title( char(plot_title(spect_type(ii),:)) );
          ylabel( 'phase' );
        endif
      endfor
    endif
  endif

end