如何在 MATLAB 中对分段周期函数进行矢量化？答案

【问题标题】：How to vectorize a piecewise periodic function in MATLAB?如何在 MATLAB 中对分段周期函数进行矢量化？
【发布时间】：2016-01-28 03:38:01
【问题描述】：

我注意到 matlab 内置函数可以处理标量或向量参数。示例：

sin(pi/2)
ans =
     1

sin([0:pi/5:pi])
ans =
         0    0.5878    0.9511    0.9511    0.5878    0.0000

如果我自己写函数，比如分段周期函数：

function v = foo(t)
t = mod( t, 2 ) ;

if ( t < 0.1 )
   v = 0 ;
elseif ( t < 0.2 )
   v = 10 * t - 1 ;
else
   v = 1 ;
end

我可以将其称为单个值：

[foo(0.1) foo(0.15) foo(0.2)]
ans =
         0    0.5000    1.0000

但是，如果函数的输入是向量，它不会像内置函数那样自动向量化：

foo([0.1:0.05:0.2])
ans =
     1

是否有可以在函数定义中使用的语法来指示如果提供了向量，则应该生成向量？或者像 sin, cos, ... 这样的内置函数会检查输入的类型，如果输入是向量，是否会产生相同的结果？

【问题讨论】：

标签： matlab

【解决方案1】：

您需要稍微更改语法才能处理任何大小的数据。我通常使用逻辑过滤器来向量化 if 语句，就像您尝试做的那样：

function v = foo(t)

v = zeros(size(t));
t = mod( t, 2 ) ;

filt1 = t<0.1;
filt2 = ~filt1 & t<0.2;
filt3 = ~filt1 & ~filt2;

v(filt1) = 0;
v(filt2) = 10*t(filt2)-1;
v(filt3) = 1;

在这段代码中，我们有三个逻辑过滤器。第一个选择所有元素，例如t<0.1。第二个挑选出所有元素，例如不在第一个过滤器中的t<0.2。最终过滤器得到其他所有内容。

然后我们使用它来设置向量v。我们将v 中与第一个过滤器匹配的每个元素设置为0。我们将v 中的所有内容都设置为与第二个过滤器匹配的10*t-1。我们将匹配第三个过滤器的v 的每个元素设置为1。

要更全面地了解矢量化，请查看上面的MATLAB help page。

【讨论】：

【解决方案2】：

一种使操作次数最少的简单方法是：

function v = foo(t)
t = mod(t, 2);
v = ones(size(t)) .* (t > 0.1);
v(t < 0.2) = 10*t(t < 0.2) - 1;
end

如果向量很大，使用ind = t < 0.2 可能会更快，并在最后一行使用它。这样，您只需在数组中搜索一次。此外，乘法可能会被带有逻辑索引的额外行替换。

【讨论】：

【解决方案3】：

我反复遇到同样的问题，因此我一直在寻找更通用的解决方案并想出了这个：

%your function definition
c={@(t)(mod(t,2))<0.1,0,...
@(t)(mod(t,2))<0.2,@(t)(10 * t - 1),...
true,1};
%call pw which returns the function
foo=pw(c{:});
%example evaluation
foo([0.1:0.05:0.2])

现在是 pw 的代码

function f=pw(varargin)
for ip=1:numel(varargin)
    switch class(varargin{ip})
        case {'double','logical'}
            varargin{ip}=@(x)(repmat(varargin{ip},size(x)));
        case 'function_handle'
            %do nothing
        otherwise
            error('wrong input class')
    end
end
c=struct('cnd',varargin(1:2:end),'fcn',varargin(2:2:end));
f=@(x)pweval(x,c);

end

function y=pweval(x,p)
todo=true(size(x));
y=x.*0;
for segment=1:numel(p)
    mask=todo;
    mask(mask)=logical(p(segment).cnd(x(mask)));
    y(mask)=p(segment).fcn(x(mask));
    todo(mask)=false;
end
assert(~any(todo));
end

【讨论】：