sgwt_cheby_square : Chebyshev coefficients for square of polynomial function d=sgwt_cheby_square(c) Inputs : c - Chebyshev coefficients for p(x) = sum c(1+k) T_k(x) ; 0<=K<=M Outputs : d - Chebyshev coefficients for p(x)^2 = sum d(1+k) T_k(x) ; 0<=k<=2*M
0001 % sgwt_cheby_square : Chebyshev coefficients for square of polynomial 0002 % 0003 % function d=sgwt_cheby_square(c) 0004 % 0005 % Inputs : 0006 % c - Chebyshev coefficients for p(x) = sum c(1+k) T_k(x) ; 0<=K<=M 0007 % 0008 % Outputs : 0009 % d - Chebyshev coefficients for p(x)^2 = sum d(1+k) T_k(x) ; 0010 % 0<=k<=2*M 0011 0012 % This file is part of the SGWT toolbox (Spectral Graph Wavelet Transform toolbox) 0013 % Copyright (C) 2010, David K. Hammond. 0014 % 0015 % The SGWT toolbox is free software: you can redistribute it and/or modify 0016 % it under the terms of the GNU General Public License as published by 0017 % the Free Software Foundation, either version 3 of the License, or 0018 % (at your option) any later version. 0019 % 0020 % The SGWT toolbox is distributed in the hope that it will be useful, 0021 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0022 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0023 % GNU General Public License for more details. 0024 % 0025 % You should have received a copy of the GNU General Public License 0026 % along with the SGWT toolbox. If not, see <http://www.gnu.org/licenses/>. 0027 0028 function d=sgwt_cheby_square(c) 0029 M=numel(c)-1; 0030 cp=c; 0031 cp(1)=.5*c(1); 0032 0033 % adjust cp so that 0034 % p(x) = sum cp(1+k) T_k(x) 0035 % for all k>=0 (rather than with special case for k=0) 0036 % 0037 % Then formula for dp in terms of cp is simpler. 0038 % Ref: my notes, july 20, 2009 0039 dp=zeros(1,2*M+1); 0040 % keep in mind : due to indexing from 1 0041 % c(1+k) is k'th Chebyshev coefficient 0042 0043 for m=0:(2*M) 0044 if (m==0) 0045 dp(1+m)=dp(1+m)+.5*cp(1)^2; 0046 for i=0:M 0047 dp(1+m)=dp(1+m)+.5*cp(1+i)^2; 0048 end 0049 elseif (m<=M) 0050 for i=0:m 0051 dp(1+m)=dp(1+m)+.5*cp(1+i)*cp(1+m-i); 0052 end 0053 for i=0:(M-m) 0054 dp(1+m)=dp(1+m)+.5*cp(1+i)*cp(1+i+m); 0055 end 0056 for i=m:M 0057 dp(1+m)=dp(1+m)+.5*cp(1+i)*cp(1+i-m); 0058 end 0059 else % M<m<=2*M 0060 for i=(m-M):M 0061 dp(1+m)=dp(1+m)+.5*cp(1+i)*cp(1+m-i); 0062 end 0063 end 0064 end 0065 d=dp; 0066 d(1)=2*dp(1);