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 Z0001 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 % See files gpl-3.0.txt and lgpl-3.0.txt included in this distribution. 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