# v_quadpeak

## PURPOSE

V_PEAK2DQUAD find quadratically-interpolated peak in a N-D array

## SYNOPSIS

function [v,x,t,m,ze]=v_quadpeak(z)

## DESCRIPTION

    V_PEAK2DQUAD find quadratically-interpolated peak in a N-D array

      Inputs:  Z(m,n,...)   is the input array (ignoring trailing singleton dimensions)
                            Note: a row vector will have 2 dimensions

     Outputs:  V        is the peak value
               X(:,1)  is the position of the peak (in fractional subscript values)
               T        is -1, 0, +1 for maximum, saddle point or minimum
               M        defines the fitted quadratic: z = [x y ... 1]*M*[x y ... 1]'
               ZE       the estimated version of Z

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- v_peak2dquad V_PEAK2DQUAD find quadratically-interpolated peak in a 2D array
- v_psycest V_PSYCEST estimate multiple psychometric functions
- v_psycestu V_PSYCESTU estimate unimodal psychometric function

## SOURCE CODE

    0001 function [v,x,t,m,ze]=v_quadpeak(z)
    0002 %V_PEAK2DQUAD find quadratically-interpolated peak in a N-D array
    0003 %
    0004 %  Inputs:  Z(m,n,...)   is the input array (ignoring trailing singleton dimensions)
    0005 %                        Note: a row vector will have 2 dimensions
    0006 %
    0007 % Outputs:  V        is the peak value
    0008 %           X(:,1)  is the position of the peak (in fractional subscript values)
    0009 %           T        is -1, 0, +1 for maximum, saddle point or minimum
    0010 %           M        defines the fitted quadratic: z = [x y ... 1]*M*[x y ... 1]'
    0011 %           ZE       the estimated version of Z
    0012
    0013 %       Copyright (C) Mike Brookes 2008
    0014 %      Version: $Id: v_quadpeak.m 10865 2018-09-21 17:22:45Z dmb $
    0015 %
    0016 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0017 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0018 %
    0019 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0020 %   This program is free software; you can redistribute it and/or modify
    0021 %   it under the terms of the GNU Lesser General Public License as published by
    0022 %   the Free Software Foundation; either version 3 of the License, or
    0023 %   (at your option) any later version.
    0024 %
    0025 %   This program is distributed in the hope that it will be useful,
    0026 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0027 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0028 %   GNU Lesser General Public License for more details.
    0029 %
    0030 %   You can obtain a copy of the GNU Lesser General Public License from
    0031 %   https://www.gnu.org/licenses/ .
    0032 %
    0033 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0034
    0035 persistent wz a
    0036 % first calculate the fixed matrix, a (can be stored if sz is constant)
    0037 sz=size(z);         % size of input array
    0038 psz=prod(sz);       % number of elements in input array
    0039 dz=numel(sz);       % number of input dimensions
    0040 mz=find(sz>1);      % non-singleton dimension indices
    0041 nm=numel(mz);       % number of non-singleton dimensions
    0042 vz=sz(mz);          % size of squeezed input array
    0043 dx=max(mz);         % number of output dimensions
    0044 if ~nm              % if the input array is a scalar
    0045     error('Cannot find peak of a scalar');
    0046 end
    0047 nc=(nm+1)*(nm+2)/2;  % number of columns in A matrix
    0048 if min(vz)<3
    0049     error('Need at least 3 points in each non-singleton dimension');
    0050 end
    0051 if isempty(wz) || numel(wz)~=numel(vz) || ~all(wz==vz)
    0052     wz=vz;
    0053     a=ones(psz,nc);
    0054     ix=(0:psz-1)';
    0055     for i=1:nm
    0056         jx=floor(ix/sz(mz(i)));
    0057         a(:,i+nc-nm-1)=1+ix-jx*sz(mz(i));
    0058         ix=jx;
    0059         a(:,(i^2-i+2)/2:i*(i+1)/2)=a(:,nc-nm:i+nc-nm-1).*repmat(a(:,i+nc-nm-1),1,i);
    0060     end
    0061     a=(a'*a)\a';        % converts to polynomial coeficients {x^2 xy y^2 x y 1]
    0062 end
    0063
    0064 % now find the peak
    0065
    0066 c=a*z(:);   % polynomial coefficients for this data
    0067 w=zeros(nm+1,nm+1);
    0068 i=1:(nm+1)*(nm+2)/2;
    0069 j=floor((sqrt(8*i-7)-1)/2);
    0070 w(i+j.*(2*nm+1-j)/2)=c;
    0071 w=(w+w.')/2; % make it symmetrical
    0072 mr=w(1:nm,1:nm);
    0073 we=w(1:nm,nm+1);
    0074 y=-(mr\we);
    0075 v=y'*we+w(nm+1,nm+1);  % value at peak
    0076
    0077 % insert singleton dimensions into outputs
    0078
    0079 x=zeros(dx,1);
    0080 x(mz)=y;
    0081 m=zeros(dx+1,dx+1);
    0082 mz(nm+1)=dx+1;
    0083 m(mz,mz)=w;
    0084 if nargout>2
    0085     ev=eig(mr);
    0086     t=all(ev>0)-all(ev<0);
    0087 end
    0088 if nargout>4
    0089     ze=zeros(sz);
    0090     scp=cumprod([1 sz(1:end-1)]);
    0091     ivec=fix(repmat((0:psz-1)',1,dz)./repmat(scp,psz,1));
    0092     xe=[1+ivec-repmat(sz,psz,1).*fix(ivec./repmat(sz,psz,1)) ones(psz,1)];
    0093     ze=reshape(sum((xe*m).*xe,2),sz);
    0094 end
    0095
    0096
    0097 if ~nargout && nm<=2
    0098     % plot the data
    0099     desc={'Maximum','Saddle Point','Minimum'};
    0100     if nargout<=2
    0101         ev=eig(mr);
    0102         t=all(ev>0)-all(ev<0);
    0103     end
    0104     if nm==1
    0105         xax=linspace(1,psz,100);
    0106         plot(xax,c(1)*xax.^2+c(2)*xax+c(3),'-r',1:psz,z(:),'ob',x,v,'^k');
    0107         set(gca,'xlim',[0.9 psz+0.1]);
    0108         ylabel('z');
    0109         xlabel(sprintf('x%d',mz(1)));
    0110         title(sprintf('\\Delta = %s: z(%.2g) = %.2g',desc{t+2},y(1),v));
    0111     else
    0112         ngr=17;
    0113         xax=repmat(linspace(1,vz(1),ngr)',1,ngr);
    0114         yax=repmat(linspace(1,vz(2),ngr),ngr,1);
    0115         zq=(c(1)*xax+c(2)*yax+c(4)).*xax+(c(3)*yax+c(5)).*yax+c(6);
    0116         hold off
    0117         mesh(xax,yax,zq,'EdgeColor','r');
    0118         hold on
    0119         plot3(repmat((1:vz(1))',1,vz(2)),repmat(1:vz(2),vz(1),1),reshape(z,vz),'ob',y(1),y(2),v,'^k');
    0120         hold off
    0121         set(gca,'xlim',[0.9 vz(1)+0.1],'ylim',[0.9 vz(2)+0.1]);
    0122         xlabel(sprintf('x%d',mz(1)));
    0123         ylabel(sprintf('x%d',mz(2)));
    0124         zlabel('z');
    0125         title(sprintf('\\Delta = %s: z(%.2g,%.2g) = %.2g',desc{t+2},y(1),y(2),v));
    0126     end
    0127 end
    0128
    0129