Home > voicebox > quadpeak.m

quadpeak

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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:

SOURCE CODE ^

0001 function [v,x,t,m,ze]=quadpeak(z)
0002 %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: quadpeak.m 713 2011-10-16 14:45:43Z 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 General Public License as published by
0022 %   the Free Software Foundation; either version 2 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 General Public License for more details.
0029 %
0030 %   You can obtain a copy of the GNU General Public License from
0031 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0032 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
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

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003