代码之家 › 专栏 › 技术社区 › Antoni Parellada

创建具有系数和可变条目数的二维余弦波矩阵

matlab

Antoni Parellada · 技术社区 · 6 年前

在我发帖后 this question yesterday ,我意识到我想用表单的每个条目创建不同n x n维度的类似矩阵。

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

其中a是系数的向量;j、x、k和y是从0到n-1的索引。

例如,如果 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

结果矩阵将有16个条目,可通过函数计算:

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)

当然,除了需要在余弦前面提供系数外,如果我想生成一个256 x 256的矩阵,那么键入所有这些余弦表达式是不可行的,例如…

我玩了for循环,但没有得到我想要的结果,在函数中获得了关于独立索引循环数的错误。

0 回复 | 直到 6 年前

HansHirse 6 年前

编辑: 我编辑了我的初始答案,添加了 Guille's comment . (一开始没有看到…)请看更新的代码。

Smee again . 您可以这样组合匿名函数/函数句柄:

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

不过,我想,有必要封装函数(句柄)的设置。 f 在一些“real”matlab函数中,请参见以下代码:

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

下面是一个测试脚本:

% 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)

输出:

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

当然,“矢量化”方法在性能方面胜出: