MATLABで学ぶ信号処理・音響工学音場制御のプログラミングについて5 （Planarity Optimization）

前回の復習
Planarityとの出会い
Planarityの理論
MATLABで論文の検証
1. 関連

前回の復習

今回紹介する音場最適化手法（Planarity）は激ムズなので、前回の記事の内容を確認してから本記事を読むことをお勧めいたいします。

MATLABで学ぶ信号処理・音響工学　音場制御のプログラミングについて1

本記事ではフィードフォワードの音場制御方法の理論について、数式とMATLABのコードを交えながら説明いたします。理論としてはすごくシンプルなので、高校生でも理解できると思います。

MATLABで学ぶ信号処理・音響工学　音場制御のプログラミングについて2 （最小二乗法を用いた最適化）

最小二乗法を用いた音場最適化制御方法について、理論式やMATLABプログラムを示しながら説明します。ここまで簡単に詳しく説明してくれている文献はほかにないので、ぜひ読んでみてください。

MATLABで学ぶ信号処理・音響工学　音場制御のプログラミングについて3 （効果的な目標音場の設定）

最小二乗法を用いた音場最適化制御方法について、理論式やMATLABプログラムを示しながら説明します。この音場最適化手法では、目標音場の設定が重要となります。この目標音場の設定方法について詳しく説明いたします。ここまで簡単に詳しく説明してくれている文献はほかにないので、ぜひ読んでみてください。

MATLABで学ぶ信号処理・音響工学　音場制御の評価指標について1

音場制御では、制御がうまくいっているか、ブライトゾーンとダークゾーンの構築がうまくいっているかを評価するための指標があります。本記事では、Array Effort、Acoustic Contrast、ブライトゾーン内の音圧のばらつき、について説明し、MATLABのプログラムを示しながら、音場制御方法の違いと評価指標の差異を確認します。

MATLABで学ぶ信号処理・音響工学　音場制御の評価指標について1

Planarityとの出会い

今回紹介するPlanarityという音場最適化手法はイギリスのサリー大学(University of Surrey)のSpeech and Signal Processing というグループが発表したものです。

Audio Engineering Societyという学会に参加したときに、この手法の発表をしていました。
2013年当時の私は「やられたな」、「0.5歩ぐらい先を行かれているな」と感じました。ただ、自慢になりますが、発表を聞いてから理論を理解し、3週間後には論文と全く同じ解析をトレースできました。（世界の先端についていけている実感がありました。学生に戻りたいな。。。。(笑)）

Planarityを簡単に説明すると、過去の記事で紹介した”Acoustic Contrast Optimization”と”Lease Squares Optimization”の良いとこ取りができそうな手法です。

論文
Philip Coleman, Philip Jackson, Marek Olik, and Abildgaard Pedersen, “Optimizing the planarity of sound zones”, AES 52nd International conference, Guildford UK, September 2-4 2013.
https://www.researchgate.net/publication/257918766_Optimizing_the_Planarity_of_Sound_Zones

ちなみに、上記の論文では”Lease Squares Optimization”のことを”Pressure matching”と呼んでいます。また、”Planarity”は”Acoustic Contrast Optimization”の派生形なので、”Acoustic Contrast with Planarity Optimization”と呼んでいる論文もあります。

Planarityの理論

先にも述べたように、理論はむずいです。なので、簡単な紹介だと思って読み進めてください。

“Planarity”はブライトゾーン内の音場が、平面波が伝搬するような音場にすることができる音場制御方法（音場最適化手法）です。

“Acoustic Contrast Optimization”はAcoustic Contrastを最大にする音場最適化手法でした。そのため、Acoustic Contrastを最大することはできますが、「音を聴く環境としてはどうなの？」という疑問があります。

ブライトゾーンで音を聴いたときに、あまりにも変な音環境（残響が大きかったり、左右の耳の音量が違ったり）することは良くありません。

この”Acoustic Contrast Optimization”の課題を平面波が伝搬するような音場を構成できるように改良しよう！というのがPlanarityのコンセプトです

“Planarity”の評価関数は下式になります。ちなみに評価関数J_ACPはスカラーです。

$J_{ACP} = q^H G_D^H G_D q + λ_{C}(q^H G_D^H H_B^H Γ H_B G_D q -B) + λ_{M}(q^H q – E_m)$

ここで、H_Bはブライトゾーンのステアリング行列、Γは平面波の角度を調節するための重み行列です。

ちなみに、”Acoustic Contrast Optimization”の評価関数J_ACは下式です。
$J_{AC} = q^H G_D^H G_D q + λ_{C}(q^H G_D^H G_D q -B) + λ_{M}(q^H q – E_m)$

“Acoustic Contrast Optimization”の評価関数にH_BとΓを足しただけですね。

ただ、この2つの行列の算出方法がむずいんすよ。ステアリング行列は下式で表されます。

$H_B=\left( \begin{array}{c} h_1 \\ \vdots \\ h_i \\ \vdots \\ h_{LB} \end{array} \right)$

ここで、hiは下記の最大固有値の固有ベクトルと一致します。

$(S_i^H S_i+βI)^{-1} P_i^H P_i$

Siはパスバンド、Piはストップバンドで、下式で表します。

$P_i=( g_{p,c} ) , S_i=( g_{s,c} )$

gi,cはグリーン関数とバンドレンジiを用いて下式で表します。

$g_{i,c} = \frac{e^{j k r_c u_i} p_B}{L_B}$
$u_i = \left( \begin{array}{c} sinφ \\ cosφ \end{array} \right)$

長くなりそうなので、とりあえず理論式はここまでにします。詳しく知りたい方はぜひ論文を読んでみてください。

MATLABで論文の検証

では、プログラムを組んで理論を検証しましょう。論文の解析条件を図1に示します。図1の○は入力（点音源）、・は応答点を示します。応答点の群で構成される右側のゾーンがブライトゾーンで、左側がダークゾーンです。

図1　論文の解析条件

解析結果を図2に示します。subplot(2,1,1-3)が音圧分布でsubplot(2,1,4-6)が位相の分布図です。この解析の目的は、ゾーンの前方から音が聞こえるように音場を制御することです。PC(Planarity Control)では前方から後方（図の上から下方向）に音が伝搬しています。ACC(Acoustic Contrast Control)では音の方向は定まっていません。PCandACC(Pressure matching Control and Acoustic Contrast Control)では似せたような解析はできますが、すこし、前方ではなく、少し斜めに音が伝搬しています。

図2　論文の検証結果

ちなみに、PCandACC(Pressure matching Control and Acoustic Contrast Control)は”Acoustic Contrast Optimization”の評価関数と重み係数を乗じた”Lease Squares Optimization”の評価関数を合体させた音場最適化手法です。

サリー大学では「PlanarityはAcoustic Contrastを高く保ちつつ、ブライトゾーンの音場環境をデザインできる手法」と主張しています。

Acoustic Contrastを高く保ててるか気になりますよね？気になる方は自分で計算してみてください。

今回はこのへんでGood luck

プログラムは下記のfunctionファイルをカレントディレクトリに保存してください。

function [G]=cell_matrix(G0,N1,N2,f)
%cell_matrix
%It transform from cell to matrix.
%for example
%G--->{2,2}(1:10)
%G{1,1}(1)=1
%G{2,1}(1)=2
%G{1,2}(1)=3
%G{2,2}(1)=4
%A=[1 3
%   2 4]
 
G=zeros(N1,N2);
 
for n1=1:N1
    for n2=1:N2
        if isnan(G0{n1,n2}(f))==1
        G(n1,n2)=0;
        else
            G(n1,n2)=G0{n1,n2}(f);
        end
         
    end
end
 
end

function  [A]=cellzeros(n1,n2,lengthnakami)
A=cell(n1,n2);
for n11=1:n1
    for n22=1:n2
        A{n11,n22}(1:lengthnakami,1)=zeros(lengthnakami,1);
    end
end
 
end

function [r]=radius(x1,x2,y1,y2,z1,z2)
 
r=sqrt( ......
    ( x2 - x1 )^2 .....
    +( y2 - y1 )^2 .....
    +( z2 - z1 )^2 .....
    );
 
end

たしか、通常の「eigs.m」ではエラーがでるので、少し書き換えているはず。

function  varargout = EIGS(varargin)
%EIGS  Find a few eigenvalues and eigenvectors of a matrix using ARPACK
%   D = EIGS(A) returns a vector of A's 6 largest magnitude eigenvalues.
%   A must be square and should be large and sparse.
%
%   [V,D] = EIGS(A) returns a diagonal matrix D of A's 6 largest magnitude
%   eigenvalues and a matrix V whose columns are the corresponding
%   eigenvectors.
%
%   [V,D,FLAG] = EIGS(A) also returns a convergence flag. If FLAG is 0 then
%   all the eigenvalues converged; otherwise not all converged.
%
%   EIGS(A,B) solves the generalized eigenvalue problem A*V == B*V*D. B
%   must be symmetric (or Hermitian) positive definite and the same size as
%   A. EIGS(A,[],...) indicates the standard eigenvalue problem A*V == V*D.
%
%   EIGS(A,K) and EIGS(A,B,K) return the K largest magnitude eigenvalues.
%
%   EIGS(A,K,SIGMA) and EIGS(A,B,K,SIGMA) return K eigenvalues. If SIGMA is:
%      'LM' or 'SM' - Largest or Smallest Magnitude
%   For real symmetric problems, SIGMA may also be:
%      'LA' or 'SA' - Largest or Smallest Algebraic
%      'BE' - Both Ends, one more from high end if K is odd
%   For nonsymmetric and complex problems, SIGMA may also be:
%      'LR' or 'SR' - Largest or Smallest Real part
%      'LI' or 'SI' - Largest or Smallest Imaginary part
%   If SIGMA is a real or complex scalar including 0, EIGS finds the
%   eigenvalues closest to SIGMA. For scalar SIGMA, and when SIGMA = 'SM',
%   B need only be symmetric (or Hermitian) positive semi-definite since it
%   is not Cholesky factored as in the other cases.
%
%   EIGS(A,K,SIGMA,OPTS) and EIGS(A,B,K,SIGMA,OPTS) specify options:
%   OPTS.issym: symmetry of A or A-SIGMA*B represented by AFUN [{false} | true]
%   OPTS.isreal: complexity of A or A-SIGMA*B represented by AFUN [false | {true}]
%   OPTS.tol: convergence: Ritz estimate residual <= tol*NORM(A) [scalar | {eps}]
%   OPTS.maxit: maximum number of iterations [integer | {300}]
%   OPTS.p: number of Lanczos vectors: K+1<p<=N [integer | {2K}]
%   OPTS.v0: starting vector [N-by-1 vector | {randomly generated}]
%   OPTS.disp: diagnostic information display level [0 | {1} | 2]
%   OPTS.cholB: B is actually its Cholesky factor CHOL(B) [{false} | true]
%   OPTS.permB: sparse B is actually CHOL(B(permB,permB)) [permB | {1:N}]
%   Use CHOL(B) instead of B when SIGMA is a string other than 'SM'.
%
%   EIGS(AFUN,N) accepts the function AFUN instead of the matrix A. AFUN is
%   a function handle and Y = AFUN(X) should return
%      A*X            if SIGMA is unspecified, or a string other than 'SM'
%      A\X            if SIGMA is 0 or 'SM'
%      (A-SIGMA*I)\X  if SIGMA is a nonzero scalar (standard problem)
%      (A-SIGMA*B)\X  if SIGMA is a nonzero scalar (generalized problem)
%   N is the size of A. The matrix A, A-SIGMA*I or A-SIGMA*B represented by
%   AFUN is assumed to be real and nonsymmetric unless specified otherwise
%   by OPTS.isreal and OPTS.issym. In all these EIGS syntaxes, EIGS(A,...)
%   may be replaced by EIGS(AFUN,N,...).
%
%   Example:
%      A = delsq(numgrid('C',15));  d1 = eigs(A,5,'SM');
%
%   Equivalently, if dnRk is the following one-line function:
%      %----------------------------%
%      function y = dnRk(x,R,k)
%      y = (delsq(numgrid(R,k))) \ x;
%      %----------------------------%
%
%      n = size(A,1);  opts.issym = 1;
%      d2 = eigs(@(x)dnRk(x,'C',15),n,5,'SM',opts);
%
%   See also EIG, SVDS, ARPACKC, FUNCTION_HANDLE.
 
%   Copyright 1984-2007 The MathWorks, Inc.
%   $Revision: 1.45.4.8 $  $Date: 2007/05/23 18:54:51 $
 
%   EIGS provides the reverse communication interface to ARPACK library
%   routines. EIGS attempts to provide an interface for as many different
%   algorithms as possible. The reverse communication interfaces are
%   documented in the ARPACK Users' Guide, ISBN 0-89871-407-9.
 
cputms = zeros(5,1);
t0 = cputime; % start timing pre-processing
       
% Process inputs and do error-checking
if (nargout > 3)
   error('MATLAB:eigs:TooManyOutputs', 'Too many output arguments.')
end
 
% Error check inputs and derive some information from them
[A,Amatrix,isrealprob,issymA,n,B,classAB,k,eigs_sigma,whch, ...
   sigma,tol,maxit,p,info,eigs_display,cholB,permB,resid,useeig,afunNargs] = ...
   checkInputs(varargin{:});
 
% Now have enough information to do early return on cases EIGS does not
% handle. For these cases, use the full EIG code.
if useeig
   fullEig(nargout);
   return
end
 
if strcmp(eigs_sigma,'SM') || ~ischar(eigs_sigma)
   % eigs(A,B,k,scalarSigma) or eigs(A,B,k,'SM'), B may be []
   % Note: sigma must be real for [s,d]saupd and [s,d]naupd
   % If sigma is complex, even if A and B are both real, we use [c,z]naupd.
   % This means that mode=3 in [s,d]naupd, which has
   % OP = real(inv(A - sigma*M)*M) and B = M
   % reduces to the same OP as [s,d]saupd and [c,z]naupd.
   % A*x = lambda*M*x, M symmetric (positive) semi-definite
   % => OP = inv(A - sigma*M)*M and B = M
   % => shift-and-invert mode
   mode = 3;
elseif isempty(B)
   % eigs(A,k,stringSigma) or eigs(A,[],k,stringSigma), stringSigma~='SM'
   % A*x = lambda*x
   % => OP = A and B = I
   mode = 1;
else
   % eigs(A,B,k,stringSigma), stringSigma~='SM'
   % A*x = lambda*B*x
   % Since we can always Cholesky factor B, follow the advice of
   % Remark 3 in ARPACK Users' Guide, and do not use mode = 2.
   % Instead, use mode = 1 with OP(x) = R'\(A*(R\x)) and B = I
   % where R is B's upper triangular Cholesky factor: B = R'*R.
   % Finally, V = R\V returns the actual generalized eigenvectors of (A,B).
   mode = 1;
end
 
if cholB || ((mode == 1) && ~isempty(B))
   % The reordering permutation permB is [] unless B is sparse
   [RB,RBT,permB] = CHOLfactorB;
end
 
permAsB = [];
if (mode == 3) && Amatrix % need lu(A-sigma*B)
   % The reordering permutation permAsB is [] unless A-sigma*B is sparse
   [L,U,P,permAsB] = LUfactorAminusSigmaB;
end % if (mode == 3) && Amatrix
 
% Allocate outputs and ARPACK work variables
if isrealprob
   if issymA % real and symmetric
      if strcmp(classAB,'single')
         aupdfun = 'ssaupd';
         eupdfun = 'sseupd';
      else
         aupdfun = 'dsaupd';
         eupdfun = 'dseupd';
      end
      lworkl = int32(p*(p+8));
      d = zeros(k,1,classAB);
   else % real but not symmetric
      if strcmp(classAB,'single')
         aupdfun = 'snaupd';
         eupdfun = 'sneupd';
      else
         aupdfun = 'dnaupd';
         eupdfun = 'dneupd';
      end
      lworkl = int32(3*p*(p+2));
      workev = zeros(3*p,1,classAB);
      d = zeros(k+1,1,classAB);
      di = zeros(k+1,1,classAB);
   end
   v = zeros(n,p,classAB);
   workd = zeros(n,3,classAB);
   workl = zeros(lworkl,1,classAB);
else % complex
   if strcmp(classAB,'single')
      aupdfun = 'cnaupd';
      eupdfun = 'cneupd';
   else
      aupdfun = 'znaupd';
      eupdfun = 'zneupd';
   end
   zv = zeros(2*n*p,1,classAB);
   workd = complex(zeros(n,3,classAB));
   zworkd = zeros(2*numel(workd),1,classAB);
   lworkl = int32(3*p^2+5*p);
   workl = zeros(2*lworkl,1,classAB);
   workev = zeros(2*2*p,1,classAB);
   zd = zeros(2*(k+1),1,classAB);
   rwork = zeros(p,1,classAB);
end
ldv = int32(n);
ipntr = zeros(15,1,'int32');
ido = int32(0); % reverse communication parameter, initial value
if isempty(B) || (mode == 1)
   bmat = 'I'; % standard eigenvalue problem
else
   bmat = 'G'; % generalized eigenvalue problem
end
nev = int32(k); % number of eigenvalues requested
ncv = int32(p); % number of Lanczos vectors
iparam = zeros(11,1,'int32');
% iparam(1) = ishift = 1 ensures we are never asked to handle ido=3
iparam([1 3 7]) = [1 maxit mode];
select = zeros(p,1,'int32');
 
% To Do: Remove this error when ARPACKC supports singles
if strcmp(classAB,'single')
   error('MATLAB:eigs:single', ...
      'EIGS does not support single precision inputs.')
end
 
% The ARPACK routines return to EIGS many times per each iteration but we
% only want to display the Ritz values once per iteration (if opts.disp>0).
% Keep track of whether we've displayed this iteration yet in eigs_iter.
eigs_iter = 0;
 
cputms(1) = cputime - t0; % end timing pre-processing
 
% Iterate until ARPACK's reverse communication parameter ido says to stop
while (ido ~= 99)
 
   t0 = cputime; % start timing ARPACK calls **aupd
 
   if isrealprob
      arpackc( aupdfun, ido, ...
         bmat, int32(n), whch, nev, tol, resid, ncv, ...
         v, ldv, iparam, ipntr, workd, workl, lworkl, info );
   else
      % The FORTRAN ARPACK routine expects the complex input zworkd to have
      % real and imaginary parts interleaved, but the OP about to be
      % applied to workd expects it in MATLAB's complex representation with
      % separate real and imaginary parts. Thus we need both.
      zworkd(1:2:end-1) = real(workd);
      zworkd(2:2:end) = imag(workd);
      arpackc( aupdfun, ido, ...
         bmat, int32(n), whch, nev, tol, resid, ncv, ...
         zv, ldv, iparam, ipntr, zworkd, workl, lworkl, rwork, info );
      workd = reshape(complex(zworkd(1:2:end-1),zworkd(2:2:end)),[n,3]);
   end
 
   if (info < 0)
      error('MATLAB:eigs:ARPACKroutineError', ...
         'Error with ARPACK routine %s: info = %d', ...
         aupdfun,full(double(info)))
   end
 
   cputms(2) = cputms(2) + (cputime-t0); % end timing ARPACK calls **aupd
   t0 = cputime; % start timing MATLAB OP(X)
 
   % Compute which columns of workd ipntr references
   cols = checkIpntr;
 
   % The ARPACK reverse communication parameter ido tells EIGS what to do
   switch ido
      case {-1,1} % abs(ido)==1 => workd(:,col2) = OP*workd(:,col1)
         switch mode
            case 1 % mode==1 => OP(x) = K*x
               if isempty(B) % standard eigenvalue problem
                  % OP(x) = A*x
                  workd(:,cols(2)) = Amtimes(workd(:,cols(1)));
               else % generalized eigenvalue problem
                  % OP(x) = R'\(A*(R\x))
                  workd(:,cols(2)) = ...
                     RBTsolve(Amtimes(RBsolve(workd(:,cols(1)))));
               end
            case 3 % mode==3 => OP(x) = inv(A-sigma*B)*B*x
               if isempty(B) % standard eigenvalue problem
                  workd(:,cols(2)) = AminusSigmaBsolve(workd(:,cols(1)));
               else % generalized eigenvalue problem
                  switch ido
                     case -1
                        workd(:,cols(2)) = Bmtimes(workd(:,cols(1)));
                        workd(:,cols(2)) = ...
                           AminusSigmaBsolve(workd(:,cols(2)));
                     case 1
                        % mode==3 and ido==1:
                        % workd(:,col2) = inv(A-sigma*B)*B*x
                        % but B*x is already pre-computed in workd(:,col3)
                        workd(:,cols(2)) = ...
                           AminusSigmaBsolve(workd(:,cols(3)));
                     otherwise
                        error('MATLAB:eigs:UnknownRCP',...
                           'Unknown reverse communication parameter.')
                  end % switch ido (inner)
               end % if isempty(B)
            otherwise % mode is not 1 or 3
               error('MATLAB:eigs:UnknownMode','Unknown mode.')
         end % switch (mode)
      case 2 % ido==2 => workd(:,col2) = B*workd(:,col1)
         if (mode == 3)
            workd(:,cols(2)) = Bmtimes(workd(:,cols(1)));
         else
            error('MATLAB:eigs:UnknownMode','Unknown mode.')
         end
      case 3 % ido==3 => EIGS does not know how to compute shifts
         % setting iparam(1) = ishift = 1 ensures this never happens
         warning('MATLAB:eigs:WorklShiftsUnsupported', ...
            ['EIGS does not support computing the shifts in workl.' ...
            ' Returning immediately.'])
         ido = int32(99);
      case 99 % ido==99 => ARPACK is done
      otherwise
         error('MATLAB:eigs:UnknownReverseCommParamFromARPACK',...
            ['Unknown value of reverse communication parameter' ...
            ' returned from %s.'],aupdfun)
 
   end % switch ido (outer)
 
   cputms(3) = cputms(3) + (cputime-t0); % end timing MATLAB OP(X)
 
   if eigs_display
%       displayRitzValues;
   end
 
end % while (ido ~= 99)
 
t0 = cputime; % start timing post-processing
 
if (info < 0)
   error('MATLAB:eigs:ARPACKroutineError', ...
      'Error with ARPACK routine %s: info = %d',aupdfun,full(info));
end % if (info < 0)
 
if (nargout >= 2)
   rvec = int32(true); % compute eigenvectors
else
   rvec = int32(false); % do not compute eigenvectors
end
 
if isrealprob
   if issymA
      arpackc( eupdfun, rvec, 'A', select, ...
         d, v, ldv, sigma, ...
         bmat, int32(n), whch, nev, tol, resid, ncv, ...
         v, ldv, iparam, ipntr, workd, workl, lworkl, info );
      if strcmp(whch,'LM') || strcmp(whch,'LA')
         d = flipud(d);
         if (rvec == 1)
            v(:,1:k) = v(:,k:-1:1);
         end
      end
      if ((strcmp(whch,'SM') || strcmp(whch,'SA')) && (rvec == 0))
         d = flipud(d);
      end
   else
      % If sigma is complex, isrealprob=true and we use [c,z]neupd.
      % So use sigmar=sigma and sigmai=0 here in dneupd.
      arpackc( eupdfun, rvec, 'A', select, ...
         d, di, v, ldv, sigma, 0, workev, ...
         bmat, int32(n), whch, nev, tol, resid, ncv, ...
         v, ldv, iparam, ipntr, workd, workl, lworkl, info );
      d = complex(d,di);
      if rvec
         d(k+1) = [];
      else
         zind = find(d == 0);
         if isempty(zind)
            d = d(k+1:-1:2);
         else
            d(max(zind)) = [];
            d = flipud(d);
         end
      end
   end
else
   zsigma = [real(sigma); imag(sigma)];
   arpackc( eupdfun, rvec, 'A', select, ...
      zd, zv, ldv, zsigma, workev, ...
      bmat, int32(n), whch, nev, tol, resid, ncv, zv, ...
      ldv, iparam, ipntr, zworkd, workl, lworkl, ...
      rwork, info );
   if issymA
      d = zd(1:2:end-1);
   else
      d = complex(zd(1:2:end-1),zd(2:2:end));
   end
   v = reshape(complex(zv(1:2:end-1),zv(2:2:end)),[n p]);
end
 
flag = processEUPDinfo(nargin<3);
 
if (issymA) || (~isrealprob)
   if (nargout <= 1)
      if isrealprob
         varargout{1} = d;
      else
         varargout{1} = d(k:-1:1,1);
      end
   else
      varargout{1} = v(:,1:k);
      varargout{2} = diag(d(1:k,1));
      if (nargout >= 3)
         varargout{3} = flag;
      end
   end
else
   if (nargout <= 1)
      varargout{1} = d;
   else
      cplxd = find(di ~= 0);
      % complex conjugate pairs of eigenvalues occur together
      cplxd = cplxd(1:2:end);
      v(:,[cplxd cplxd+1]) = [complex(v(:,cplxd),v(:,cplxd+1)) ...
         complex(v(:,cplxd),-v(:,cplxd+1))];
      varargout{1} = v(:,1:k);
      varargout{2} = diag(d);
      if (nargout >= 3)
         varargout{3} = flag;
      end
   end
end
 
if (nargout >= 2) && (mode == 1) && ~isempty(B)
   varargout{1} = RBsolve(varargout{1});
end
 
cputms(4) = cputime-t0; % end timing post-processing
 
cputms(5) = sum(cputms(1:4)); % total time
 
if (eigs_display == 2)
   printTimings;
end
 
%-------------------------------------------------------------------------%
% Nested functions
%-------------------------------------------------------------------------%
 
% checkInputs error checks the inputs to EIGS and also derives some
%   variables from them:
% A may be a matrix or a function applying OP.
% Amatrix is true if A is a matrix, false if A is a function.
% isrealprob is true if all of A, B and sigma are real, false otherwise.
% issymA is true if A is symmetric, false otherwise.
% n is the size of (square) A and B.
% B is [] for the standard problem. Otherwise it may be one of B, CHOL(B)
%   or CHOL(B(permB,permB)).
% classAB is single if either A or B is single, otherwise double.
% k is the number of eigenvalues to be computed.
% eigs_sigma is the value for sigma passed in by the user, 'LM' if it was
%   unspecified. eigs_sigma may be either a string or a scalar value.
% whch is the ARPACK string corresponding to eigs_sigma and mode.
% sigma is the ARPACK scalar corresponding to eigs_sigma and mode.
% tol is the convergence tolerance.
% maxit is the maximum number of iterations.
% p is the number of Lanczos vectors.
% info is the start value, initialized to 1 or 0 to indicate whether to use
% resid as the start vector or not.
% eigs_display is true if Ritz values should be displayed, false otherwise.
% cholB is true if CHOL(B) was passed in instead of B, false otherwise.
% permB may be [], otherwise it is the permutation in CHOL(B(permB,permB)).
% resid is the start vector if specified and info=1, otherwise all zero.
% useeig is true if we need to use EIG instead of ARPACK, otherwise false.
% afunNargs is the range of EIGS' varargin that are to be passed as
%   trailing parameters to the function as in afun(X,P1,P2,...).
   function [A,Amatrix,isrealprob,issymA,n,B,classAB,k, ...
         eigs_sigma,whch,sigma,tol,maxit,p,info,eigs_display,cholB,...
         permB,resid,useeig,afunNargs] = checkInputs(varargin)
      % Process inputs and do error-checking
 
      % Process the input A or the inputs AFUN and N
      % Start to derive some qualities (real, symmetric) about the problem
      if isfloat(varargin{1})
         A = varargin{1};
         Amatrix = true;
      else
         % By checking the function A with fcnchk, we can now use direct
         % function evaluation on the result, without resorting to feval
         A = fcnchk(varargin{1});
         Amatrix = false;
      end
      % isrealprob = isreal(A) && isreal(B) && isreal(sigma)
      isrealprob = true;
      issymA = false;
      if Amatrix
         isrealprob = isreal(A);
         issymA = ishermitian(A);
         [m,n] = size(A);
         if (m ~= n)
            error('MATLAB:eigs:NonSquareMatrixOrFunction',...
               'A must be a square matrix or a function.')
         end
      else
         n = varargin{2};
         nstr = 'Size of problem, ''n'', must be a positive integer.';
         if ~isscalar(n) || ~isreal(n)
            error('MATLAB:eigs:NonPosIntSize', nstr)
         end
         if issparse(n)
            n = full(n);
         end
         if (round(n) ~= n)
            warning('MATLAB:eigs:NonPosIntSize',['%s\n         ' ...
               'Rounding input size.'],nstr)
            n = round(n);
         end
      end
 
      % Process the input B and derive the class of the problem.
      % Is B present in the eigs call or not?
      Bpresent = true;
      Bstr = ['Generalized matrix B must be the same size as A and' ...
         ' either a symmetric positive (semi-)definite matrix or' ...
         ' its Cholesky factor.'];
      if (nargin < (3-Amatrix))
         B = [];
         Bpresent = false;
      else
         % Is the next input B or K?
         B = varargin{3-Amatrix};
         if ~isempty(B) % allow eigs(A,[],k,sigma,opts);
            if isscalar(B)
               if n ~= 1
                  % this input is really K and B is not specified
                  B = [];
                  Bpresent = false;
               else
                  % This input could be B or K.
                  % If A is scalar, then the only valid value for k is 1.
                  % So if this input is scalar, let it be B, namely
                  % eigs(4,2,...) assumes A=4, B=2, NOT A=4, k=2
                  if ~isnumeric(B)
                     error('MATLAB:eigs:BsizeMismatchAorNotSPDorNotChol', Bstr);
                  end
                  % Unless, of course, the scalar is 1, in which case 
                  % assume the that it is meant to be K.
                  if (B == 1) && ((Amatrix && nargin <= 3) || ...
                         (~Amatrix && nargin <= 4))
                      B = [];
                      Bpresent = false;
                  elseif ~isfloat(B)
                     error('MATLAB:eigs:BsizeMismatchAorNotSPDorNotChol', Bstr);
                  end
               end
            else
               % B is a not a scalar.
               if ~isfloat(B) || ~isequal(size(B),[n,n])
                  error('MATLAB:eigs:BsizeMismatchAorNotSPDorNotChol', Bstr);
               end
               isrealprob = isrealprob && isreal(B);
            end
         end
      end
      % ARPACK can only handle homogeneous inputs
      if Amatrix
         classAB = superiorfloat(A,B);
         A = cast(A,classAB);
         B = cast(B,classAB);
      else
         if ~isempty(B)
            classAB = class(B);
         else
            classAB = 'double';
         end
      end
       
      % argOffset tells us where to get the eigs inputs K, SIGMA and OPTS.
      % If A is really the function afun, then it also helps us find the
      % trailing parameters in eigs(afun,n,[B],k,sigma,opts,P1,P2,...)
      % Values of argOffset:
      %  0: Amatrix is false and Bpresent is true:
      %     eigs(afun,n,B,k,sigma,opts,P1,P2,...)
      %  1: Amatrix and Bpresent are both true, or both false
      %     eigs(A,B,k,sigma,opts)
      %     eigs(afun,n,k,sigma,opts,P1,P2,...)
      %  2: Amatrix is true and Bpresent is false:
      %     eigs(A,k,sigma,opts)
      argOffset = Amatrix + ~Bpresent;
 
      if Amatrix && ((nargin - Bpresent)>4)
         error('MATLAB:eigs:TooManyInputs', 'Too many inputs.')
      end
 
      % Process the input K.
      if (nargin < (4-argOffset))
         k = min(n,6);
      else
         k = varargin{4-argOffset};
         kstr = ['Number of eigenvalues requested, k, must be a' ...
            ' positive integer <= n.'];
         if ~isnumeric(k) || ~isscalar(k) || ~isreal(k) || (k>n)
            error('MATLAB:eigs:NonIntegerEigQty', kstr)
         end
         if issparse(k)
            k = full(k);
         end
         if (round(k) ~= k)
            warning('MATLAB:eigs:NonIntegerEigQty',['%s\n         ' ...
               'Rounding number of eigenvalues.'],kstr)
            k = round(k);
         end
      end
 
      % Process the input SIGMA and derive ARPACK values whch and sigma.
      % eigs_sigma is the value documented in the help as "SIGMA" that is
      % passed in to EIGS. eigs_sigma may be either a scalar, including 0,
      % or a string, including 'SM'.
      % In ARPACK, eigs_sigma corresponds to two variables:
      % 1.  which, called "whch" to avoid conflict with MATLAB's function
      % 2.  sigma
      % whch is always a string. sigma is always a scalar.
      % Valid combinations are shown below. Note eigs_sigma = 0/'SM' has
      % the same sigma/whch values as eigs_sigma='LM' (default) so these
      % must be distinguished by the mode.
      % eigs_sigma = 'SM' or 0 => sigma = 0, whch = 'LM' (mode=3)
      % eigs_sigma is a string not 'SM' => sigma = 0, whch = eigs_sigma (mode=1)
      % eigs_sigma is any scalar => sigma = eigs_sigma, whch = 'LM'
      % (mode=1)
      whchstr = 'Eigenvalue range sigma must be a valid 2-element string.';
      if (nargin < (5-argOffset))
         % default: eigs 'LM' => ARPACK which='LM', sigma=0
         eigs_sigma = 'LM';
         whch = 'LM';
         sigma = 0;
      else
         eigs_sigma = varargin{5-argOffset};
         if ischar(eigs_sigma)
            % eigs(string) => ARPACK which=string, sigma=0
            if ~isequal(size(eigs_sigma),[1,2])
               error('MATLAB:eigs:EigenvalueRangeNotValid', ...
                  [whchstr '\nFor real symmetric A, the' ...
                  ' choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.' ...
                  '\nFor non-symmetric or complex' ...
                  ' A, the choices are ''%s'', ''%s'', ''%s'', ''%s'',' ...
                  ' ''%s'' or ''%s''.\n'], ...
                  'LM','SM','LA','SA','BE','LM','SM','LR','SR','LI','SI')
            end
            eigs_sigma = upper(eigs_sigma);
            if strcmp(eigs_sigma,'SM')
               % eigs('SM') => ARPACK which='LM', sigma=0
               whch = 'LM';
            else
               % eigs(string), where string~='SM' => ARPACK which=string, sigma=0
               whch = eigs_sigma;
            end
            sigma = zeros(classAB);
         else
            % eigs(scalar) => ARPACK which='LM', sigma=scalar
            if ~isfloat(eigs_sigma) || ~isscalar(eigs_sigma)
               error('MATLAB:eigs:EigenvalueShiftNonScalar',...
                  'Eigenvalue shift sigma must be a scalar.')
            end
            sigma = eigs_sigma;
            if issparse(sigma)
               sigma = full(sigma);
            end
            sigma = cast(sigma,classAB);
            isrealprob = isrealprob && isreal(sigma);
            whch = 'LM';
         end
      end
 
      % Process the input OPTS and derive some ARPACK values.
      % ARPACK's minimum tolerance is eps/2 ([S/D]LAMCH's EPS)
      tol = eps(classAB);
      maxit = [];
      p = [];
      % Always use resid as the start vector, whether it is OPTS.v0 or
      % randomly generated within eigs.  We default resid to empty here.
      % If the user does not initialize it, we provide a random residual
      % below.
      info = int32(1);
      resid = []; 
      eigs_display = 1;
      cholB = false; % do we have B or its Cholesky factor?
      permB = []; % if cholB, is it chol(B), or chol(B(permB,permB))?
      if (nargin >= (6-argOffset))
         opts = varargin{6-argOffset};
         if ~isa(opts,'struct')
            error('MATLAB:eigs:OptionsNotStructure',...
               'Options argument must be a structure.')
         end
         if isfield(opts,'issym') && ~Amatrix
            issymA = opts.issym;
            if (issymA ~= false) && (issymA ~= true)
               error('MATLAB:eigs:InvalidOptsIssym', ...
                  'opts.issym must be true or false.')
            end
         end
         if isfield(opts,'isreal') && ~Amatrix
            if (opts.isreal ~= false) && (opts.isreal ~= true)
               error('MATLAB:eigs:InvalidOptsIsreal', ...
                  'opts.isreal must be true or false.')
            end
            isrealprob = isrealprob && opts.isreal;
         end
         if ~isempty(B) && (isfield(opts,'cholB') || isfield(opts,'permB'))
            if isfield(opts,'cholB')
               cholB = opts.cholB;
               if (cholB ~= false) && (cholB ~= true)
                  error('MATLAB:eigs:InvalidOptsCholB', ...
                     'opts.cholB must be true or false.')
               end
               if isfield(opts,'permB')
                  if issparse(B) && cholB
                     permB = opts.permB;
                     if ~isvector(permB) || ~isequal(sort(permB(:)),(1:n)')
                        error('MATLAB:eigs:InvalidOptsPermB',...
                           'opts.permB must be a permutation of 1:n.')
                     end
                  else
                     warning('MATLAB:eigs:IgnoredOptionPermB', ...
                        ['Ignoring opts.permB since B is not its sparse' ...
                        ' Cholesky factor.'])
                  end
               end
            end
         end
         if isfield(opts,'tol')
            if ~isfloat(tol) || ~isscalar(opts.tol) || ~isreal(opts.tol) || (opts.tol<=0)
               error('MATLAB:eigs:InvalidOptsTol',...
                  ['Convergence tolerance opts.tol must be a strictly' ...
                  ' positive real scalar.'])
            end
            tol = cast(full(opts.tol),classAB);
         end
         if isfield(opts,'p')
            p = opts.p;
            pstr = ['Number of basis vectors opts.p must be a positive' ...
               ' integer <= n.'];
            if ~isnumeric(p) || ~isscalar(p) || ~isreal(p) || (p<=0) || (p>n)
               error('MATLAB:eigs:InvalidOptsP', pstr)
            end
            if issparse(p)
               p = full(p);
            end
            if (round(p) ~= p)
               warning('MATLAB:eigs:NonIntegerVecQty',['%s\n         ' ...
                  'Rounding number of basis vectors.'],pstr)
               p = round(p);
            end
         end
         if isfield(opts,'maxit')
            maxit = opts.maxit;
            str = ['Maximum number of iterations opts.maxit must be' ...
               ' a positive integer.'];
            if ~isnumeric(maxit) || ~isscalar(maxit) || ~isreal(maxit) || (maxit<=0)
               error('MATLAB:eigs:OptsMaxitNotPosInt', str)
            end
            if issparse(maxit)
               maxit = full(maxit);
            end
            if (round(maxit) ~= maxit)
               warning('MATLAB:eigs:NonIntegerIterationQty',['%s\n         ' ...
                  'Rounding number of iterations.'],str)
               maxit = round(maxit);
            end
         end
         if isfield(opts,'v0')
            if ~isfloat(opts.v0) || ~isequal(size(opts.v0),[n,1])
               error('MATLAB:eigs:WrongSizeOptsV0',...
                  'Start vector opts.v0 must be n-by-1.')
            end
            if isrealprob
               if ~isreal(opts.v0)
                  error('MATLAB:eigs:NotRealOptsV0',...
                     'Start vector opts.v0 must be real for real problems.')
               end
               resid(1:n,1) = full(opts.v0);
            else
               resid(2:2:2*n,1) = full(imag(opts.v0));
               resid(1:2:(2*n-1),1) = full(real(opts.v0));
            end
         end
         if isfield(opts,'disp')
            eigs_display = opts.disp;
            dispstr = 'Diagnostic level opts.disp must be an integer.';
            if ~isnumeric(eigs_display) || ~isscalar(eigs_display) || ...
                  ~isreal(eigs_display) || (eigs_display<0)
               error('MATLAB:eigs:NonIntegerDiagnosticLevel', dispstr)
            end
            if (round(eigs_display) ~= eigs_display)
               warning('MATLAB:eigs:NonIntegerDiagnosticLevel', ...
                  '%s\n         Rounding diagnostic level.',dispstr)
               eigs_display = round(eigs_display);
            end
         end
         if isfield(opts,'cheb')
            error('MATLAB:eigs:ObsoleteOptionCheb', ...
               'Polynomial acceleration opts.cheb is an obsolete option.');
         end
         if isfield(opts,'stagtol')
            error('MATLAB:eigs:ObsoleteOptionStagtol', ...
               'Stagnation tolerance opts.stagtol is an obsolete option.');
         end
      end
      if (isempty(resid))
        if isrealprob
           resid = cast(rand(n,1),classAB);
        else
           resid = cast(rand(2*n,1),classAB);
        end
      end
 
      afunNargs = zeros(1,0);
      if ~Amatrix
         % The trailing parameters for afun start at varargin{7-argOffset}
         % in eigs(afun,n,[B],k,sigma,opts,P1,P2,...). If there are no
         % trailing parameters in eigs, then afunNargs is a 1-by-0 empty
         % and no trailing parameters are passed to afun(x)
         afunNargs = 7-argOffset:nargin;
      end
 
      % Now that OPTS has been processed, do final error checking and
      % assign ARPACK variables
 
      % Extra check on input B
      if ~isempty(B)
         % B must be symmetric (Hermitian) positive (semi-)definite
         if cholB
            if ~isequal(triu(B),B)
               error('MATLAB:eigs:BsizeMismatchAorNotSPDorNotChol', Bstr)
            end
         else
            if ~ishermitian(B)
               error('MATLAB:eigs:BsizeMismatchAorNotSPDorNotChol', Bstr)
            end
         end
      end
 
      % Extra check on input K
      % We fall back on using the full EIG code if K is too large.
      useeig = false;
      if isrealprob && issymA
         knstr = sprintf(['For real symmetric problems, must have' ...
            ' number of eigenvalues k < n.\n']);
      else
         knstr = sprintf(['For nonsymmetric and complex problems,' ...
            ' must have number of eigenvalues k < n-1.\n']);
      end
      if isempty(B)
         knstr = [knstr 'Using eig(full(A)) instead.'];
      else
         knstr = [knstr 'Using eig(full(A),full(B)) instead.'];
      end
      if (k == 0)
         useeig = true;
      end
      if isrealprob && issymA
         if (k > n-1)
            if (n >= 6)
               warning('MATLAB:eigs:TooManyRequestedEigsForRealSym', ...
                  '%s',knstr)
            end
            useeig = true;
         end
      else
         if (k > n-2)
            if (n >= 7)
               warning('MATLAB:eigs:TooManyRequestedEigsForComplexNonsym', ...
                  '%s',knstr)
            end
            useeig = true;
         end
      end
 
      % Extra check on input SIGMA
      if isrealprob && issymA
         if ~isreal(sigma)
            error('MATLAB:eigs:ComplexShiftForRealSymProblem',...
               ['For real symmetric problems, eigenvalue shift sigma must' ...
               ' be real.'])
         end
      else
         if ~isrealprob && issymA && ~isreal(sigma)
            warning('MATLAB:eigs:ComplexShiftForHermitianProblem', ...
               ['Complex eigenvalue shift sigma on a Hermitian problem' ...
               ' (all real eigenvalues).'])
         end
      end
      if isrealprob && issymA
         if strcmp(whch,'LR')
            whch = 'LA';
            warning('MATLAB:eigs:SigmaChangedToLA', ...
               ['For real symmetric problems, sigma value ''LR''' ...
               ' (Largest Real) is now ''LA'' (Largest Algebraic).'])
         end
         if strcmp(whch,'SR')
            whch = 'SA';
            warning('MATLAB:eigs:SigmaChangedToSA', ...
               ['For real symmetric problems, sigma value ''SR''' ...
               ' (Smallest Real) is now ''SA'' (Smallest Algebraic).'])
         end
         if ~ismember(whch,{'LM', 'SM', 'LA', 'SA', 'BE'})
            error('MATLAB:eigs:EigenvalueRangeNotValid', ...
               [whchstr '\nFor real symmetric A, the' ...
               ' choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ...
               'LM','SM','LA','SA','BE');
         end
      else
         if strcmp(whch,'BE')
            warning('MATLAB:eigs:SigmaChangedToLM', ...
               ['Sigma value ''BE'' is now only available for real' ...
               ' symmetric problems.  Computing ''LM'' eigenvalues instead.'])
            whch = 'LM';
         end
         if ~ismember(whch,{'LM', 'SM', 'LR', 'SR', 'LI', 'SI'})
            error('MATLAB:eigs:EigenvalueRangeNotValid', ...
               [whchstr '\nFor non-symmetric or complex' ...
               ' A, the choices are ''%s'', ''%s'', ''%s'', ''%s'',' ...
               ' ''%s'' or ''%s''.\n'],'LM','SM','LR','SR','LI','SI');
         end
      end
       
      % The remainder of the error checking does not apply for the large
      % values of K that force us to use full EIG instead of ARPACK.
      if useeig
         return
      end
 
      % Extra check on input OPTS.p
      if isempty(p)
         if isrealprob && ~issymA
            p = min(max(2*k+1,20),n);
         else
            p = min(max(2*k,20),n);
         end
      else
         if isrealprob && issymA
            if (p <= k)
               error('MATLAB:eigs:InvalidOptsPforRealSymProb',...
                  ['For real symmetric problems, must have number of' ...
                  ' basis vectors opts.p > k.'])
            end
         else
            if (p <= k+1)
               error('MATLAB:eigs:InvalidOptsPforComplexOrNonSymProb',...
                  ['For nonsymmetric and complex problems, must have number of' ...
                  ' basis vectors opts.p > k+1.'])
            end
         end
      end
 
      % Extra check on input OPTS.maxit
      if isempty(maxit)
         maxit = max(300,ceil(2*n/max(p,1)));
      end
 
   end % checkInputs
 
%-------------------------------------------------------------------------%
   function fullEig(nOutputs)
      % Use EIG(FULL(A)) or EIG(FULL(A),FULL(B)) instead of ARPACK
      if ~isempty(B)
         B = Bmtimes(eye(n));
      end
      if isfloat(A)
         if issparse(A);
            A = full(A);
         end
      else
         % A is specified by a function.
         % Form the matrix A by applying the function.
         if ischar(eigs_sigma) && ~strcmp(eigs_sigma,'SM')
            % A is a function multiplying A*x
            AA = eye(n);
            for i = 1:n
               AA(:,i) = A(AA(:,i),varargin{afunNargs});
            end
            A = AA;
         else
            if (isfloat(eigs_sigma) && eigs_sigma == 0) || strcmp(eigs_sigma,'SM')
               % A is a function solving A\x
               invA = eye(n);
               for i = 1:n
                  invA(:,i) = A(invA(:,i),varargin{afunNargs});
               end
               A = eye(n) / invA;
            else
               % A is a function solving (A-sigma*B)\x
               % B may be [], indicating the identity matrix
               % U = (A-sigma*B)\sigma*B
               % => (A-sigma*B)*U = sigma*B
               % => A*U = sigma*B(U + eye(n))
               % => A = sigma*B(U + eye(n)) / U
               if isempty(B)
                  sB = eigs_sigma*eye(n);
               else
                  sB = eigs_sigma*B;
               end
               U = zeros(n,n);
               for i = 1:n
                  U(:,i) = A(sB(:,i),varargin{afunNargs});
               end
               A = sB*(U+eye(n)) / U;
            end
         end
      end
 
      if isempty(B)
         eigInputs = {A};
      else
         eigInputs = {A,B};
      end
      % Now with full floating point matrices A and B, use EIG:
      if (nOutputs <= 1)
         d = eig(eigInputs{:});
      else
         [V,D] = eig(eigInputs{:});
         d = diag(D);
      end
 
      % Grab the eigenvalues we want, based on sigma
      firstKindices = 1:k;
      lastKindices = n:-1:n-k+1;
      if ischar(eigs_sigma)
         switch eigs_sigma
            case 'LM'
               [ignore,ind] = sort(abs(d));
               range = lastKindices;
            case 'SM'
               [ignore,ind] = sort(abs(d));
               range = firstKindices;
            case 'LA'
               [ignore,ind] = sort(d);
               range = lastKindices;
            case 'SA'
               [ignore,ind] = sort(d);
               range = firstKindices;
            case 'LR'
               [ignore,ind] = sort(abs(real(d)));
               range = lastKindices;
            case 'SR'
               [ignore,ind] = sort(abs(real(d)));
               range = firstKindices;
            case 'LI'
               [ignore,ind] = sort(abs(imag(d)));
               range = lastKindices;
            case 'SI'
               [ignore,ind] = sort(abs(imag(d)));
               range = firstKindices;
            case 'BE'
               [ignore,ind] = sort(abs(d));
               range = [1:floor(k/2), n-ceil(k/2)+1:n];
            otherwise
               error('MATLAB:eigs:fullEigSigma','Unknown value of sigma');
         end
      else
         % sigma is a scalar
         [ignore,ind] = sort(abs(d-eigs_sigma));
         range = 1:k;
      end
       
      if (nOutputs <= 1)
         varargout{1} = d(ind(range));
      else
         varargout{1} = V(:,ind(range));
         varargout{2} = D(ind(range),ind(range));
         if (nOutputs == 3)
            % flag indicates "convergence"
            varargout{3} = 0;
         end
      end
       
   end % FULLEIG
    
%-------------------------------------------------------------------------%
   function [RB,RBT,perm] = CHOLfactorB
      % permB may be [] (from checkInputs) if the problem is not sparse
      % or if it was not passed in as opts.permB
      perm = permB;
      if cholB
         % CHOL(B) was passed in as B
         RB = B;
         RBT = B';
      else
         % CHOL(B) was not passed into EIGS
         if (mode == 1) && ~isempty(B)
            % Algorithm requires CHOL(B) to be computed
            if issparse(B)
               perm = symamd(B);
               [RB,pB] = chol(B(perm,perm));
            else
               [RB,pB] = chol(B);
            end
            if (pB == 0)
               RBT = RB';
            else
               error('MATLAB:eigs:BNotSPD', ...
                  'B is not symmetric positive definite.')
            end
         end
      end
   end % CHOLfactorB
 
%-------------------------------------------------------------------------%
   function [L,U,P,perm] = LUfactorAminusSigmaB
      % LU factor A-sigma*B, including a reordering perm if it is sparse
      if isempty(B)
         if issparse(A)
            AsB = A - sigma * speye(n);
         else
            AsB = A - sigma * eye(n);
         end
      else
         if cholB
            if issparse(B)
               AsB = A - sigma * Bmtimes(speye(n));
            else
               AsB = A - sigma * Bmtimes(eye(n));
            end
         else
            AsB = A - sigma * B;
         end
      end
      if issparse(AsB)
         [L,U,P,Q] = lu(AsB);
         [perm,ignore] = find(Q);
      else
         [L,U,P] = lu(AsB);
         perm = [];
      end
      % Warn if lu(A-sigma*B) is ill-conditioned
      % => sigma is close to an exact eigenvalue of (A,B)
      dU = diag(U);
      rcondestU = full(min(abs(dU)) / max(abs(dU)));
      if (rcondestU < eps)
         if isempty(B)
            ds = '(A-sigma*I)';
         else
            ds = '(A-sigma*B)';
         end
         warning('MATLAB:eigs:SigmaNearExactEig',...
            [ds ' has small reciprocal condition' ...
            ' estimate: %f\n' ...
            '         indicating that sigma is near an exact' ...
            ' eigenvalue.\n         The algorithm may not converge unless' ...
            ' you try a new value for sigma.\n'], ...
            rcondestU);
      end
   end % LUfactorAminusSigmaB
 
%-------------------------------------------------------------------------%
   function cols = checkIpntr
      % Check that ipntr returned from ARPACK refers to the start of a
      % column of workd.
      if ~isempty(B) && (mode == 3) && (ido == 1)
         inds = double(ipntr(1:3));
      else
         inds = double(ipntr(1:2));
      end
      [rows,cols] = ind2sub([n,3],inds);
      nonOneRows = find(rows~=1);
      if ~isempty(nonOneRows)
         error('MATLAB:eigs:ipntrMismatchWorkdColumn', ...
         ['One of ipntr(1:3) does not refer to the start' ...
            ' of a column of the %d-by-3 array workd.'],n)
      end
   end % checkIpntr
 
%-------------------------------------------------------------------------%
   function v = Amtimes(u)
      % Matrix-vector multiply v = A*u
      if Amatrix
         v = A * u;
      else % A is a function
         v = A(u,varargin{afunNargs});
         if isrealprob && ~isreal(v)
            error('MATLAB:eigs:complexFunction', ...
                  'AFUN is complex; set opts.isreal = false.');
         end
      end
   end
 
%-------------------------------------------------------------------------%
   function v = Bmtimes(u)
      % Matrix-vector multiply v = B*u
      if cholB % use B's cholesky factor and its transpose
         if ~isempty(permB)
            v(permB,:) = RBT * (RB * u(permB,:));
         else
            v = RBT * (RB * u);
         end
      else
         v = B * u;
      end
   end
 
%-------------------------------------------------------------------------%
   function v = RBsolve(u)
      % Solve v = RB\u for v
      if issparse(B)
         if ~isempty(permB)
            v(permB,:) = RB \ u;
         else
            v = RB \ u;
         end
      else
         RBopts.UT = true;
         v = linsolve(RB,u,RBopts);
      end
   end
 
%-------------------------------------------------------------------------%
   function v = RBTsolve(u)
      % Solve v = RB'\u for v
      if issparse(B)
         if ~isempty(permB)
            v = RBT \ u(permB,:);
         else
            v = RBT \ u;
         end
      else
         RBTopts.LT = true;
         v = linsolve(RBT,u,RBTopts);
      end
   end
 
%-------------------------------------------------------------------------%
   function v = AminusSigmaBsolve(u)
      % Solve v = (A-sigma*B)\u for v
      if Amatrix
         if ~isempty(permAsB)
            % use LU reordering permAsB
            v(permAsB,:) = U \ (L \ (P * u));
         else
            v = U \ (L \ (P * u));
         end
      else % A is a function
         v = A(u,varargin{afunNargs});
         if isrealprob && ~isreal(v)
            error('MATLAB:eigs:complexFunction', ...
                  'AFUN is complex; set opts.isreal = false.');
         end
      end
   end % AminusSigmaBsolve
 
%-------------------------------------------------------------------------%
%    function displayRitzValues
%       % Display a few Ritz values at the current iteration
%       iter = double(ipntr(15));
%       if (iter > eigs_iter) && (ido ~= 99)
%          eigs_iter = iter;
% %          ds = sprintf(['Iteration %d: a few Ritz values of the' ...
% %             ' %d-by-%d matrix:'],iter,p,p);
% %          disp(ds)
%          if isrealprob
%             if issymA
%                dispvec = workl(double(ipntr(6))+(0:p-1));
%                if strcmp(whch,'BE')
%                   % roughly k Large eigenvalues and k Small eigenvalues
%                   disp(dispvec(max(end-2*k+1,1):end))
%                else
%                   % k eigenvalues
%                   disp(dispvec(max(end-k+1,1):end))
%                end
%             else
%                dispvec = complex(workl(double(ipntr(6))+(0:p-1)), ...
%                   workl(double(ipntr(7))+(0:p-1)));
%                % k+1 eigenvalues (keep complex conjugate pairs together)
%                disp(dispvec(max(end-k,1):end))
%             end
%          else
%             dispvec = complex(workl(2*double(ipntr(6))-1+(0:2:2*(p-1))), ...
%                workl(2*double(ipntr(6))+(0:2:2*(p-1))));
%             disp(dispvec(max(end-k+1,1):end))
%          end
%       end
%    end
 
%-------------------------------------------------------------------------%
   function flag = processEUPDinfo(warnNonConvergence)
      % Process the info flag returned by the ARPACK routine **eupd
      flag = 0;
      if (info ~= 0)
         es = ['Error with ARPACK routine ' eupdfun ':\n'];
         switch double(info)
            case 2
               ss = sum(select);
               if (ss < k)
                  error('MATLAB:eigs:ARPACKroutineError02ssLTk', ...
                     [es 'The logical variable select was only set' ...
                     ' with %d 1''s instead of nconv=%d (k=%d).\n' ...
                     'Please report this to the ARPACK authors at' ...
                     ' arpack@caam.rice.edu.'], ...
                     ss,double(iparam(5)),k)
               else
                  error('MATLAB:eigs:ARPACKroutineError02', ...
                     [es 'The LAPACK reordering routine %strsen' ...
                     ' did not return all %d eigenvalues.'], ...
                     aupdfun(1),k);
               end
            case 1
               error('MATLAB:eigs:ARPACKroutineError01', ...
                  [es 'The Schur form could not be reordered by the' ...
                  ' LAPACK routine %strsen.\nPlease report this to the' ...
                  ' ARPACK authors at arpack@caam.rice.edu.'], ...
                  aupdfun(1))
            case -14
               error('MATLAB:eigs:ARPACKroutineErrorMinus14', ...
                  [es aupdfun ...
                  ' did not find any eigenvalues to sufficient accuracy.']);
            otherwise
               error('MATLAB:eigs:ARPACKroutineError', ...
                  [es 'info = %d. Please consult the ARPACK Users''' ...
                  ' Guide for more information.'],full(info));
         end
      else
         nconv = double(iparam(5));
         if (nconv == 0)
            if (warnNonConvergence)
               warning('MATLAB:eigs:NoEigsConverged', ...
                  'None of the %d requested eigenvalues converged.',k)
            else
               flag = 1;
            end
         elseif (nconv < k)
            if (warnNonConvergence)
               warning('MATLAB:eigs:NotAllEigsConverged', ...
                  'Only %d of the %d requested eigenvalues converged.', ...
                  nconv,k)
            else
               flag = 1;
            end
         end
      end
   end % processEUPDinfo
 
%-------------------------------------------------------------------------%
   function printTimings
      % Print the time taken for each major stage of the EIGS algorithm
      if (mode == 1)
         innerstr = sprintf(['Compute A*X:' ...
            '                               %f\n'],cputms(3));
      elseif (mode == 3)
         if isempty(B)
            innerstr = sprintf(['Solve (A-SIGMA*I)*X=Y for X:' ...
               '               %f\n'],cputms(3));
         else
            innerstr = sprintf(['Solve (A-SIGMA*B)*X=B*Y for X:' ...
               '             %f\n'],cputms(3));
         end
      end
      if ((mode == 3) && (Amatrix))
         if isempty(B)
            prepstr = sprintf(['Pre-processing, including lu(A-sigma*I):' ...
               '   %f\n'],cputms(1));
         else
            prepstr = sprintf(['Pre-processing, including lu(A-sigma*B):' ...
               '   %f\n'],cputms(1));
         end
      else
         prepstr = sprintf(['Pre-processing:' ...
            '                            %f\n'],cputms(1));
      end
      sstr = sprintf('***********CPU Timing Results in seconds***********');
      ds = sprintf(['\n' sstr '\n' ...
         prepstr ...
         'ARPACK''s %s:                           %f\n' ...
         innerstr ...
         'Post-processing with ARPACK''s %s:      %f\n' ...
         '***************************************************\n' ...
         'Total:                                     %f\n' ...
         sstr '\n'], ...
         aupdfun,cputms(2),eupdfun,cputms(4),cputms(5));
      disp(ds)
   end % printTimings
 
%-------------------------------------------------------------------------%
% End of nested functions
%-------------------------------------------------------------------------%
 
end % EIGS
 
%-------------------------------------------------------------------------%
% Subfunctions
%-------------------------------------------------------------------------%
function tf = ishermitian(A)
%ISHERMITIAN
tf = isequal(A,A');
end % ishermititan
%-------------------------------------------------------------------------%
% End of subfunctions
%-------------------------------------------------------------------------%

function [Hi]=steering_matrix(ii,beta,Stop,Pass,ZoneA,CenterZoneA,k)
 
if Stop+1<=ii&&ii<=360-Stop
    Stop_theta=[1:ii-(Stop+1) ii+(Stop+1):360];
elseif ii<Stop+1
    Stop_theta=ii+(Stop+1):1:360-(Stop-ii+1);
else
    %ii>360-Stop
    Stop_theta=1+Stop-(360-ii):1:ii-(Stop+1);    
end
 
 
if Pass+1<=ii&&ii<=360-Pass
    Pass_theta=ii-Pass:ii+Pass;
elseif ii<Pass+1
    Pass_theta=[1:ii+Pass 360-(Pass-ii):360];
else
    %ii>360-Pass
    Pass_theta=[1:Pass-(360-ii) ii-Pass:360];
end
 
D=zeros(360-(2*Stop+1),size(ZoneA,1));
B=zeros(2*Pass+1,size(ZoneA,1));
 
 
for nM=1:size(ZoneA)
    for nst=1:length(Stop_theta)
%         ust=[sin(Stop_theta(nst)/180*pi)
%             cos(Stop_theta(nst)/180*pi)];
        ust=[cos(Stop_theta(nst)/180*pi)
            sin(Stop_theta(nst)/180*pi)
            ];
%         D(nst,nM)=exp(j*k*(ZoneA(nM,:) )*ust)/size(ZoneA,1);
%         D(nst,nM)=exp(j*k*(ZoneA(nM,:)-CenterZoneA)*ust)/size(ZoneA,1);
        D(nst,nM)=exp(j*k*(ZoneA(nM,:))*ust)/size(ZoneA,1);
    end
     
    for npa=1:length(Pass_theta)
        upa=[cos(Pass_theta(npa)/180*pi)
            sin(Pass_theta(npa)/180*pi)
            ];
%         B(npa,nM)=exp(j*k*(ZoneA(nM,:) )*upa)/size(ZoneA,1);
%         B(npa,nM)=exp(j*k*(ZoneA(nM,:)-CenterZoneA)*upa)/size(ZoneA,1);
        B(npa,nM)=exp(j*k*(ZoneA(nM,:))*upa)/size(ZoneA,1);
    end
end
 
[h d]=EIGS(inv( D'*D+beta*eye(size(D'*D)) )*B'*B,1,'lm');
Hi=h.';
 
end

ここからが実行ファイル
実行ファイル1：伝達関数行列の計算

clear all
 
% f_max=10^4;
rho=1.205;
c=343;
freq=100:100:7000;
% freq=1000;
x=-0.1260:0.021:0.1260;
y=-0.1155:0.021:0.1155;
 
xy=zeros(length(x)*length(y),2);
n=0;
for X=1:length(x)
    for Y=1:length(y)
        n=1+n;
        xy(n,1:3)=[x(X) y(Y) 0];
    end
end
 
ZoneA=xy+[-0.6*ones(size(xy,1),1) zeros(size(xy,1),1) zeros(size(xy,1),1)];
ZoneB=xy+[+0.6*ones(size(xy,1),1) zeros(size(xy,1),1) zeros(size(xy,1),1)];
 
xi=[ZoneA;ZoneB];
% % % % y_secondary_i=(ys1,ys2,ys3)
R=1.2;
theta=linspace(0,2*pi-2*pi/48,48);
Speaker=[R*cos(theta).' R*sin(theta).' zeros(size(theta,2),1)];
 
Z=cellzeros(size(xi,1),size(Speaker,1),length(freq));
% Z=cellzeros(size(xi,1),size(Speaker,1),f_max);
h = waitbar(0,'Please wait...');
for f=1:1:length(freq)
% for f=1:f_max
%     w=2*pi*f;
%     k=2*pi*f/c;
    w=2*pi*freq(f);
    k=2*pi*freq(f)/c;
    for ns=1:length(Speaker)
        for nm=1:size(xi,1)
          [r]=radius(Speaker(ns,1),xi(nm,1),Speaker(ns,2),xi(nm,2),Speaker(ns,3),xi(nm,3));
          Z{nm,ns}(f)=1i*k*c*rho/(4*pi*r)*exp(-1i*k*r);
        end
    end
       waitbar(f/length(freq),h)
%        waitbar(f/f_max,h)
end
close(h)
 
save('Z100_7000.mat','Z','-v7.3')

実行ファイル2：ステアリング行列の計算

% % % % 26/09/2013
% % % % Optimizing the planarity of sound zones
clear all
 
x=-0.1260:0.021:0.1260;
y=-0.1155:0.021:0.1155;
% freq=1000;
freq=100:100:7000;
xy=zeros(length(x)*length(y),2);
n=0;
for X=1:length(x)
    for Y=1:length(y)
        n=1+n;
        xy(n,1:2)=[x(X) y(Y)];
    end
end
 
CenterZoneA=[-0.6 0];
ZoneA=xy+[-0.6*ones(size(xy,1),1) zeros(size(xy,1),1)];
ZoneB=xy+[+0.6*ones(size(xy,1),1) zeros(size(xy,1),1)];
xi=[ZoneA;ZoneB];
% % % % y_secondary_i=(ys1,ys2,ys3)
R=1.2;
theta=linspace(0,2*pi-2*pi/48,48);
 
Speaker=[R*cos(theta).' R*sin(theta).'];
 
% % % % % % % % % % % % % 
c=343;
beta=10^-4;
Stop=6;
Pass=3;
% Phi=90;
% % % % % % % % % % % % % 
 
H=cellzeros(size(Speaker,1),size(ZoneA,1),length(freq));
h = waitbar(0,'Please wait...');
for f=1:1:length(freq)
    w=2*pi*freq(f);
    k=2*pi*freq(f)/c;
    for ii=1:360
        [Hi]=steering_matrix(ii,beta,Stop,Pass,ZoneA,CenterZoneA,k);
        for n=1:size(ZoneA,1)
            H{ii,n}(f)=Hi(n).';
        end
    end
       waitbar(f/length(freq),h)
end
close(h)
save('H100_7000.mat','H','-v7.3')
 
% [q d]=eigs(inv( Gd'*Gd+beta*eye(size(D'*D)) )*(Gb'*H'*Gamma*H*Gb),1,'lm');

実行ファイル3：Planarityの入力qの計算

% % % % 30/09/2013
% % % % Optimizing the planarity of sound zones
clear all
 
load('H100_7000.mat')
load('Z100_7000.mat')
lambda2=0;
% lambda2=10^10;
x=-0.1260:0.021:0.1260;
y=-0.1155:0.021:0.1155;
freq=100:100:7000;
xy=zeros(length(x)*length(y),2);
n=0;
for X=1:length(x)
    for Y=1:length(y)
        n=1+n;
        xy(n,1:2)=[x(X) y(Y)];
    end
end
CenterZoneA=[-0.6 0];
ZoneA=xy+[-0.6*ones(size(xy,1),1) zeros(size(xy,1),1)];
ZoneB=xy+[+0.6*ones(size(xy,1),1) zeros(size(xy,1),1)];
 
xi=[ZoneA;ZoneB];
% % % % y_secondary_i=(ys1,ys2,ys3)
R=1.2;
theta=linspace(0,2*pi-2*pi/48,48);
Speaker=[R*cos(theta).' R*sin(theta).'];
 
% % % % % % % % % % % % % 
c=343;
beta=10^-4;
% % % % % % weighting
gamma=zeros(360,360);
AAA=30:1:150;
% AAA=90;
for n=1:length(AAA)
    gamma(AAA(n),AAA(n))=1;
end
% AAA=60:1:120;
% BBB=hamming(length(AAA));
% for n=1:length(AAA)
% gamma(AAA(n),AAA(n))=BBB(n);
% end
% % % % % % % % % % % % % 
 
 
C=zeros(length(freq),1);
AE=zeros(length(freq),1);
q_optimal=zeros(size(Speaker,1),length(freq));
h = waitbar(0,'Please wait...');
 
for f=1:1:length(freq)
    [G]=cell_matrix(Z,size(xi,1),size(Speaker,1),f);
    [Hb]=cell_matrix(H,360,size(ZoneA,1),f);    
    Gb=G(1:size(xi,1)/2,:);
    Gd=G(size(xi,1)/2+1:end,:);
    [q d]=eigs(inv( Gd'*Gd+lambda2*eye(size(Gd'*Gd)) )*(Gb'*Hb'*gamma*Hb*Gb),1,'lm');
    q_optimal(1:size(Speaker,1),f)=q;
    ref=mean(abs(Gb*q))/mean(abs(Gb*ones(size(Gb,2),1)));
    AE(f,1)=(q'*q)/(size(Speaker,1)*ref^2);
    C(f,1)=size(Gd,1)*q'*Gb'*Gb*q/........
        ( size(Gb,1)*q'*Gd'*Gd*q );
       waitbar(f/length(freq),h)
end
close(h)
 
save('PC100_7000.mat','AE','C','q_optimal')

実行ファイル4：表示

clear all
 
load('PC1000_optpart3.mat')
qpc=q_optimal;
% load('ACC1000_optpart3.mat')
% qacc=q_optimal;
% load('PMandACC1000_optpart3.mat')
% qpmandacc=q_optimal;
 
rho=1.205;
c=343;
R=1.2;
freq=1000;
% % % % % % 
x=-0.1260:0.021:0.1260;
y=-0.1155:0.021:0.1155;
xy=zeros(length(x)*length(y),2);
n=0;
for X=1:length(x)
    for Y=1:length(y)
        n=1+n;
        xy(n,1:2)=[x(X) y(Y)];
    end
end
CenterZoneA=[-0.6 0];
ZoneA=xy+[-0.6*ones(size(xy,1),1) zeros(size(xy,1),1)];
ZoneB=xy+[+0.6*ones(size(xy,1),1) zeros(size(xy,1),1)];
clear x
clear y
% % % % % % 
theta=linspace(0,2*pi-2*pi/48,48);
Speaker=[R*cos(theta).' R*sin(theta).'];
% % % % % % 
 
% % % % % % % % % % % % % % % % % % % % % % % % % % 
x=-1.8:0.02:1.8;
y=-1.8:0.02:1.8;
 
 
% Pressure=zeros(length(x),length(y));
% for f=1:length(freq)
for f=1
    pcPressure=zeros(length(y),length(x));
%     accPressure=zeros(length(y),length(x));
%     pmandaccPressure=zeros(length(y),length(x));
    k=2*pi*freq(f)/c;
for ns=1:size(Speaker,1)
    for X=1:length(x)
        for Y=1:length(y)
            r=sqrt( (x(X)-Speaker(ns,1))^2+(y(Y)-Speaker(ns,2))^2 );
            if r<=0.05
                r=0.05;
            end
            pcPressure(Y,X)=pcPressure(Y,X)+qpc(ns,f)/max(abs(qpc))*1i*k*c*rho/(4*pi*r)*exp(-1i*k*r);
            accPressure(Y,X)=accPressure(Y,X)+qacc(ns,f)/max(abs(qacc))*1i*k*c*rho/(4*pi*r)*exp(-1i*k*r);
            pmandaccPressure(Y,X)=pmandaccPressure(Y,X)+qpmandacc(ns,f)/max(abs(qpmandacc))*1i*k*c*rho/(4*pi*r)*exp(-1i*k*r);
 
        end
    end
end
figure(3)
subplot(231)
surf(x,y,10*log10(abs(pcPressure)))
shading interp
view([0 90])
title('PC')
colorbar
hold on
plot3(Speaker(:,1),Speaker(:,2),1000*ones(size(Speaker,1)),'wo')
plot3(ZoneA(:,1),ZoneA(:,2),10^30*ones(size(ZoneA,1)),'w.')
plot3(ZoneB(:,1),ZoneB(:,2),10^30*ones(size(ZoneB,1)),'w.')
hold off
subplot(234)
surf(x,y,angle(pcPressure))
% colormap gray
shading interp
view([0 90])
title(freq(f))
colorbar
hold on
plot3(Speaker(:,1),Speaker(:,2),1000*ones(size(Speaker,1)),'wo')
plot3(ZoneA(:,1),ZoneA(:,2),10^30*ones(size(ZoneA,1)),'w.')
plot3(ZoneB(:,1),ZoneB(:,2),10^30*ones(size(ZoneB,1)),'w.')
hold off
 
% subplot(232)
% surf(x,y,10*log10(abs(accPressure)))
% shading interp
% view([0 90])
% title('ACC')
% colorbar
% hold on
% plot3(Speaker(:,1),Speaker(:,2),1000*ones(size(Speaker,1)),'wo')
% plot3(ZoneA(:,1),ZoneA(:,2),10^30*ones(size(ZoneA,1)),'w.')
% plot3(ZoneB(:,1),ZoneB(:,2),10^30*ones(size(ZoneB,1)),'w.')
% hold off
% subplot(235)
% surf(x,y,angle(accPressure))
% % colormap gray
% shading interp
% view([0 90])
% colorbar
% title(freq(f))
% hold on
% plot3(Speaker(:,1),Speaker(:,2),1000*ones(size(Speaker,1)),'wo')
% plot3(ZoneA(:,1),ZoneA(:,2),10^30*ones(size(ZoneA,1)),'w.')
% plot3(ZoneB(:,1),ZoneB(:,2),10^30*ones(size(ZoneB,1)),'w.')
% hold off
% 
% subplot(233)
% surf(x,y,10*log10(abs(pmandaccPressure)))
% shading interp
% view([0 90])
% title('PMandACC')
% colorbar
% hold on
% plot3(Speaker(:,1),Speaker(:,2),1000*ones(size(Speaker,1)),'wo')
% plot3(ZoneA(:,1),ZoneA(:,2),10^30*ones(size(ZoneA,1)),'w.')
% plot3(ZoneB(:,1),ZoneB(:,2),10^30*ones(size(ZoneB,1)),'w.')
% hold off
% subplot(236)
% surf(x,y,angle(pmandaccPressure))
% % colormap gray
% shading interp
% view([0 90])
% title('PMandACC')
% colorbar
% hold on
% plot3(Speaker(:,1),Speaker(:,2),1000*ones(size(Speaker,1)),'wo')
% plot3(ZoneA(:,1),ZoneA(:,2),10^30*ones(size(ZoneA,1)),'w.')
% plot3(ZoneB(:,1),ZoneB(:,2),10^30*ones(size(ZoneB,1)),'w.')
% hold off
 
 
end