home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / GENPHASE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.6 KB  |  82 lines
  1. % function resp = genphase(d)
  2. %
  3. %   Uses the complex-cepstrum (Oppenheim & Schafer, Digital
  4. %   Signal Processing, pg 501) to generate a complex 
  5. %   frequency response, RESP, whose magnitude is equal to
  6. %   the real, positive response D, but whose phase corresponds
  7. %   to a stable, minimum phase function. Both D and RESP are
  8. %   VARYING matrices.
  9. %
  10. %   See also: FITMAG, FITSYS, MAGFIT, MUSYNFIT, and MUSYNFLP.
  11.  
  12. function resp = genphase(d)
  13.   if nargin ~= 1
  14.     disp(['usage: resp = genphase(d)']);
  15.     return
  16.   end
  17.   [dtype,drows,dcols,dnum] = minfo(d);
  18.   if dtype ~= 'vary' | drows ~= 1 | dcols ~= 1
  19.     error('data should be a 1x1, varying matrix')
  20.     return
  21.   end
  22.   magorig = d(1:dnum,1).';          %ROW vector
  23.   omegaorig = d(1:dnum,2).';        %ROW vector
  24.   if any(diff(omegaorig)<=0)
  25.     error('frequency data should be monotonically increasing')
  26.     return
  27.   end
  28.   pw = sqrt(omegaorig(dnum)*omegaorig(1));
  29.      [come,mag] = flatten(omegaorig,magorig);
  30.   capt = 1/pw;
  31.   tcomes = capt*capt*(come .^2);            % top half disk
  32. %
  33. % warp the frequencies around the disk
  34. % domeorig has values between 0 and pi
  35. %
  36.   domeorig = acos( (1 - tcomes) ./ ( 1 + tcomes) ); % top half disk
  37. %
  38.   lowf=domeorig(1);
  39.   highf = domeorig(length(domeorig));
  40.   if lowf <= 0 | highf >= pi
  41.     disp('frequency range should be on positive imaginary')
  42.     return
  43.   end
  44.  
  45.   smgap = min([lowf pi-highf]);
  46.   nn = ceil( (log(2*pi/smgap))/log(2) );
  47.   npts = 2*(2^nn);
  48.   hnpts = 2^nn;
  49.   if npts > 17000
  50.     disp(['large frequency range - calculation may be slow!'])
  51.   end
  52.  
  53. % interpolate the frequency and magnitude data to
  54. %    hnpts points linearly spaced around top half of disk
  55.  
  56.   lindomee = (pi/hnpts)*(0:hnpts-1);
  57.   [lindome,linmag] = terpol(domeorig,mag,hnpts);
  58. % duplicate data around disk
  59.   dome = [lindome (2*pi)-fliplr(lindome)];          %all disk
  60.   mag = [linmag  fliplr(linmag)];                   %all disk
  61.  
  62. % complex cepstrum to get min-phase
  63.  
  64.   ymag = log( mag .^2 );
  65.   ycc = ifft(ymag);                              % 2^N
  66.   nptso2 = npts/2;                               % 2^(N-1)
  67.   xcc = ycc(1:nptso2);                           % 2^(N-1)
  68.   xcc(1) = xcc(1)/2;                             % halve it at 0
  69.   xhat = exp(fft(xcc));                          % 2^(N-1)
  70.   domeg = dome(1:2:nptso2-1);                    % 2^(N-2)
  71.   xhat = xhat(1:nptso2/2);                       % 2^(N-2)
  72. % interpolate back to original frequency data
  73.   [compd,err] = terpolb(domeg,xhat,domeorig);
  74.   if err < 0
  75.     error('not sampled high enough')
  76.     return
  77.   end
  78.  
  79.   resp = vpck(compd.',omegaorig.');
  80. %
  81. % Copyright MUSYN INC 1991,  All Rights Reserved
  82.