toroidal_spect

 [amp, theta, phi] = toroidal_spect(data, dims);

 Computes coefficients using least-squares of terms for spherical harmonics of
 a function on a toroidal surface. 

 f(theta,phi) = acc(m,n)*cos(m*theta)*cos(n*phi) +
                acs(m,n)*cos(m*theta)*sin(n*phi) +
                asc(m,n)*sin(m*theta)*cos(n*phi) +
                ass(m,n)*sin(m*theta)*sin(n*phi) 
         
  NOTE: Each term is a summation over 1:m,1:n, where m and n are the order of
        the highest spherical harmonic along each dimension.

  data: the surface function, arranged as a row vector with dimensions of
  theta and phi given in dims. If data is a matrix, rows are individual
  surfaces that are to be modelled.

 
 03/25/02 Petr Janata

0001 function [amp, theta, phi, names, Rsqr, X] = toroidal_spect(data, dims,m, n)
0002 % [amp, theta, phi] = toroidal_spect(data, dims);
0003 %
0004 % Computes coefficients using least-squares of terms for spherical harmonics of
0005 % a function on a toroidal surface.
0006 %
0007 % f(theta,phi) = acc(m,n)*cos(m*theta)*cos(n*phi) +
0008 %                acs(m,n)*cos(m*theta)*sin(n*phi) +
0009 %                asc(m,n)*sin(m*theta)*cos(n*phi) +
0010 %                ass(m,n)*sin(m*theta)*sin(n*phi)
0011 %
0012 %  NOTE: Each term is a summation over 1:m,1:n, where m and n are the order of
0013 %        the highest spherical harmonic along each dimension.
0014 %
0015 %  data: the surface function, arranged as a row vector with dimensions of
0016 %  theta and phi given in dims. If data is a matrix, rows are individual
0017 %  surfaces that are to be modelled.
0018 %
0019 %
0020 % 03/25/02 Petr Janata
0021 
0022 nunits_theta = dims(1);
0023 nunits_phi = dims(2);
0024 nunits_total = nunits_theta * nunits_phi;
0025 
0026 if any(size(data) ==  1)
0027   data = data(:)';
0028   nsamps = 1;
0029   if length(data) ~= nunits_total
0030     error('toroidal_spect: length of data does not match product of dimensions')
0031   end
0032 else
0033   nsamps = size(data,1);
0034   if size(data,2) ~= nunits_total
0035     error('toroidal_spect: data size dimension mismatch')
0036   end
0037 end
0038 
0039 
0040 %
0041 % Compute each of the sin and cos components of the expansion
0042 %
0043 clear cos_theta sin_theta cos_phi sin_phi
0044 
0045 theta = 2*pi*(0:nunits_theta-1)/nunits_theta;
0046 phi = 2*pi*(0:nunits_phi-1)/nunits_phi;
0047 
0048 for iharm = 1:max([m n])
0049   harm = iharm-1;
0050   tmp = ...
0051       repmat(cos(harm*theta)',1,nunits_phi);
0052   cos_theta(:,iharm) = tmp(:);
0053   
0054   tmp = ...
0055       repmat(sin(harm*theta)',1,nunits_phi);
0056   sin_theta(:,iharm) = tmp(:);
0057   
0058   tmp = repmat(cos(harm*phi),nunits_theta,1);
0059   cos_phi(:,iharm) = tmp(:);
0060   
0061   tmp = repmat(sin(harm*phi),nunits_theta,1);
0062   sin_phi(:,iharm) = tmp(:);
0063 end
0064 
0065 %
0066 % Compute the terms of the expansion for which amplitude parameters will be estimated
0067 %
0068 clear coscos cossin sincos sinsin
0069 for mharm = 1:m
0070   for nharm = 1:n
0071     coscos(mharm,nharm,:) = cos_theta(:,mharm) .* cos_phi(:,nharm);
0072     coscos_names{mharm,nharm} = sprintf('cc%d%d',mharm-1,nharm-1);
0073       
0074     cossin(mharm,nharm,:) = cos_theta(:,mharm) .* sin_phi(:,nharm);
0075     cossin_names{mharm,nharm} = sprintf('cs%d%d',mharm-1,nharm-1);
0076 
0077     sincos(mharm,nharm,:) = sin_theta(:,mharm) .* cos_phi(:,nharm);
0078     sincos_names{mharm,nharm} = sprintf('sc%d%d',mharm-1,nharm-1);
0079 
0080     sinsin(mharm,nharm,:) = sin_theta(:,mharm) .* sin_phi(:,nharm);
0081     sinsin_names{mharm,nharm} = sprintf('ss%d%d',mharm-1,nharm-1);
0082   end % for jharm=
0083 end % for iharm=
0084 
0085 % Reshape the matrices so it is easy to access the whole set of harmonics for
0086 % any given output unit on the map.
0087 coscos = reshape(coscos,m*n,nunits_total)';
0088 cossin = reshape(cossin,m*n,nunits_total)';
0089 sincos = reshape(sincos,m*n,nunits_total)';
0090 sinsin = reshape(sinsin,m*n,nunits_total)';
0091 
0092 % Save names containing indices of each of the terms
0093 names = [coscos_names(:)' cossin_names(:)' sincos_names(:)' sinsin_names(:)'];
0094 
0095 % Put them together into a single regression matrix
0096 X = [coscos cossin sincos sinsin];
0097 
0098 % Compute some costly matrices right now
0099 pinvX = pinv(X'*X);
0100 pinvXX = pinvX*X';
0101 
0102 amp = zeros(nsamps,m*n*4);
0103 
0104 % Estimate regression coefficients for each of the samples
0105 for isamp = 1:nsamps
0106   y = data(isamp,:)';
0107   y = y-mean(y);
0108   b = pinvXX*y;
0109   y_prime = X*b;
0110   y_err = y - y_prime;
0111   Rsqr(isamp) = 1- (y_err'*y_err)/(y'*y);
0112   amp(isamp,:) = b';
0113 end