% function resp = genphase(d)
%
%   Uses the complex-cepstrum (Oppenheim & Schafer, Digital
%   Signal Processing, pg 501) to generate a complex 
%   frequency response, RESP, whose magnitude is equal to
%   the real, positive response D, but whose phase corresponds
%   to a stable, minimum phase function. Both D and RESP are
%   VARYING matrices.
%
%   See also: FITMAG, FITSYS, MAGFIT, MUSYNFIT, and MUSYNFLP.

function resp = genphase(d)
  if nargin ~= 1
    disp(['usage: resp = genphase(d)']);
    return
  end
  [dtype,drows,dcols,dnum] = minfo(d);
  if dtype ~= 'vary' | drows ~= 1 | dcols ~= 1
    error('data should be a 1x1, varying matrix')
    return
  end
  magorig = d(1:dnum,1).';          %ROW vector
  omegaorig = d(1:dnum,2).';        %ROW vector
  if any(diff(omegaorig)<=0)
    error('frequency data should be monotonically increasing')
    return
  end
  pw = sqrt(omegaorig(dnum)*omegaorig(1));
	 [come,mag] = flatten(omegaorig,magorig);
  capt = 1/pw;
  tcomes = capt*capt*(come .^2);            % top half disk
%
% warp the frequencies around the disk
% domeorig has values between 0 and pi
%
  domeorig = acos( (1 - tcomes) ./ ( 1 + tcomes) ); % top half disk
%
  lowf=domeorig(1);
  highf = domeorig(length(domeorig));
  if lowf <= 0 | highf >= pi
    disp('frequency range should be on positive imaginary')
    return
  end

  smgap = min([lowf pi-highf]);
  nn = ceil( (log(2*pi/smgap))/log(2) );
  npts = 2*(2^nn);
  hnpts = 2^nn;
  if npts > 17000
    disp(['large frequency range - calculation may be slow!'])
  end

% interpolate the frequency and magnitude data to
%    hnpts points linearly spaced around top half of disk

  lindomee = (pi/hnpts)*(0:hnpts-1);
  [lindome,linmag] = terpol(domeorig,mag,hnpts);
% duplicate data around disk
  dome = [lindome (2*pi)-fliplr(lindome)];          %all disk
  mag = [linmag  fliplr(linmag)];                   %all disk

% complex cepstrum to get min-phase

  ymag = log( mag .^2 );
  ycc = ifft(ymag);                              % 2^N
  nptso2 = npts/2;                               % 2^(N-1)
  xcc = ycc(1:nptso2);                           % 2^(N-1)
  xcc(1) = xcc(1)/2;                             % halve it at 0
  xhat = exp(fft(xcc));                          % 2^(N-1)
  domeg = dome(1:2:nptso2-1);                    % 2^(N-2)
  xhat = xhat(1:nptso2/2);                       % 2^(N-2)
% interpolate back to original frequency data
  [compd,err] = terpolb(domeg,xhat,domeorig);
  if err < 0
    error('not sampled high enough')
    return
  end

  resp = vpck(compd.',omegaorig.');
%
% Copyright MUSYN INC 1991,  All Rights Reserved
