Creating a matrix of 2D cosines waves with coefficients and variable number of entries

Question

After I posted this question yesterday, I realized that I want to create similar matrices of different n x n dimensions with each entry of the form

a * cos(j * x + k * y)

where a is a vector of coefficients; and j, x, k and y are indexes from 0 to n - 1.

If, for instance, n = 4,

>> n = 4;
>> x = 0:(n-1);
>> y = 0:(n-1);
>> [x,y] = meshgrid(x,y)
x =

   0   1   2   3
   0   1   2   3
   0   1   2   3
   0   1   2   3

y =

   0   0   0   0
   1   1   1   1
   2   2   2   2
   3   3   3   3

The resultant matrix would have 16 entries which could be computed by the function:

f = @(x, y,a0,a1,a2,a3,b0,b1,b2,b3,c0,c1,c2,c3,d0,d1,d2,d3)... 
a0*cos(0*x + 0*y) + a1*cos(0*x + 1*y) +...
a2*cos(0*x + 2*y) + a3*cos(0*x + 3*y) + ...
b0*cos(1*x + 0*y) + b1*cos(1*x + 1*y) + ...
b2*cos(1*x + 2*y) + b3*cos(1*x + 3*y) + ...
c0*cos(2*x + 1*y) + c1*cos(2*x + 1*y) + ...
c2*cos(2*x + 2*y) + c3*cos(2*x + 3*y) + ...
d0*cos(3*x + 1*y) + d1*cos(3*x + 1*y) + ...
d2*cos(3*x + 2*y) + d3*cos(3*x + 3*y)

Of course, aside from the need to furnish the coefficients in front of the cosines, typing all these cosine expressions is not doable if I want to generate a 256 x 256 matrix, for example...

I played with for-loops but I didn't get what I am after, getting error regarding the number of independent indexing loops within a function.

Answer 1

EDIT: I edited my initial answer, adding the idea given in Guille's comment. (Haven't seen that in first place...) Please, see the updated code.

---

Smee again. You can combine anonymous functions / function handles like this:

f = @(x) sin(x);
g = @(x) cos(x);
h = @(x) f(x) + g(x);

Nevertheless, I guess, it's necessary to encapsulate the setup of your function (handle) f into some "real" MATLAB function, see the following code:

function f = setupF(n, a)

  % Possibly, add some checks, e.g. for numel(a) == n^2, and so on.

  % Initialize function handle.  
  f = @(x, y) 0;
  ind = 0;

  % Iteratively add cosine parts. 
  for ii = 0:(n-1)
    for jj = 0:(n-1)
      ind = ind + 1;
      g = @(x, y) a(ind) * cos(ii * x + jj * y);
      f = @(x, y) f(x, y) + g(x, y); 
    end
  end

end

Here comes a test script:

% Set up parameters.
n = 3;
a = reshape(1:n^2, n, n);

% Set up f(x, y) by function.
f = setupF(n, a);

% Set up f explicitly, as g(x, y). 
g = @(x, y) ...
  a(1) * cos(0*x + 0*y) + ...
  a(2) * cos(0*x + 1*y) + ...
  a(3) * cos(0*x + 2*y) + ...
  a(4) * cos(1*x + 0*y) + ...
  a(5) * cos(1*x + 1*y) + ...
  a(6) * cos(1*x + 2*y) + ...
  a(7) * cos(2*x + 0*y) + ...
  a(8) * cos(2*x + 1*y) + ...
  a(9) * cos(2*x + 2*y);

% Set up f(x, y) by vectorization, as h(x, y).
I = 0:(n-1);
J = 0:(n-1);
[I, J] = meshgrid(I, J);
h = @(x, y, n, a) sum(reshape(a .* cos(x * I + y * J), n^2, 1));
h = @(x, y, n, a) arrayfun(@(x, y) h(x, y, n, a), x, y);

% Set up test data.
x = linspace(0, 2*pi, 5);
y = linspace(0, 2*pi, 5);
[X, Y] = meshgrid(x, y);

% Compare outputs.
fRet = f(X, Y)
gRet = g(X, Y)
hRet = h(X, Y, n, a)

And, the output:

fRet =
   45.0000  -18.0000   15.0000  -18.0000   45.0000
   -6.0000   -5.0000   -2.0000    5.0000   -6.0000
   15.0000   -6.0000    5.0000   -6.0000   15.0000
   -6.0000    5.0000   -2.0000   -5.0000   -6.0000
   45.0000  -18.0000   15.0000  -18.0000   45.0000

gRet =
   45.0000  -18.0000   15.0000  -18.0000   45.0000
   -6.0000   -5.0000   -2.0000    5.0000   -6.0000
   15.0000   -6.0000    5.0000   -6.0000   15.0000
   -6.0000    5.0000   -2.0000   -5.0000   -6.0000
   45.0000  -18.0000   15.0000  -18.0000   45.0000

hRet =
   45.0000  -18.0000   15.0000  -18.0000   45.0000
   -6.0000   -5.0000   -2.0000    5.0000   -6.0000
   15.0000   -6.0000    5.0000   -6.0000   15.0000
   -6.0000    5.0000   -2.0000   -5.0000   -6.0000
   45.0000  -18.0000   15.0000  -18.0000   45.0000

And, of course, the "vectorization" approach wins in terms of performance: