我正在使用简单的阈值来检测信号中的零边缘。由于噪声可能非常强烈,因此我正在使用滞后来改善结果。尽管这有很大帮助,但它也会减慢很多。有没有办法改善这一点,甚至是一种更好的方法来计算这一点?目前,我正在使用直线循环方法:
% generate a signal
t = linspace(0, 10, 1000);
y = sin(2 * pi * t);
% constants
threshold = 0;
hyst = 0.1;
% thresholding
yth = zeros(size(y));
for i = 2:length(y)
if yth(i - 1) > 0.5
yth(i) = y(i) > (threshold - hyst);
else
yth(i) = y(i) > (threshold + hyst);
end
end
与yth = y > threshold
相比,这要慢得多。
可以通过预先计算yth(i)
的两种可能性来获得改进(仅适用于大型输入;在R2017B上(:
function yth = idea1()
yth = false(size(y)); % change #1 - `false` vs `zeros`
c1 = y > (th - hyst); % change #2 - precompute c1, c2
c2 = y > (th + hyst);
for k = 2:numel(y)
if yth(k - 1)
yth(k) = c1(k);
else
yth(k) = c2(k);
end
end
end
通常,可以通过矢量化获得很大的改进。要介绍此问题,我们必须了解切换逻辑,因此让我们将其插入文字:
- 从开始模式2 。
- 模式1 :从
c1
获取值。只要c1
包含true
值 - 拿走。当遇到false
值时,切换到模式2 。 - 模式2 :从
c2
获取值。只要c2
包含false
值 - 取之量即可。当遇到true
值时,切换到模式1 。
因此,如果我们能找到过渡地点,我们实际上就完成了。
经过一定的反复试验,虽然我无法以提高性能的方式摆脱循环,但我确实得出了可能的结论(idea2
(。此外,查看编码的正确yth
,我想出了一些相当不错的近似值(idea3
和idea4
( - 尽管这些需要对其他输入进行调整。
也许有人可以使用它来找到一个更聪明的实现,并具有较少的冗余计算。我的完整代码在下面提供。编码算法的RLE是根据此答案进行了调整的,并从这里解码了。
function q48637952
% generate a signal
t = linspace(0, 10, 1000).';
y = sin(2 * pi * t);
% constants
th = 0;
hyst = 0.1;
%% Comaprison:
% Correctness:
R = {originalIdea(), idea1(), idea2()};
assert(isequal(R{:}));
% Runtime:
T = [timeit(@originalIdea,1), timeit(@idea1,1), timeit(@idea2,1)];
disp(T);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function yth = originalIdea()
yth = zeros(size(y));
for i = 2:length(y)
if yth(i - 1) > 0.5
yth(i) = y(i) > (th - hyst);
else
yth(i) = y(i) > (th + hyst);
end
end
end
function yth = idea1()
yth = false(size(y)); % change #1 - `false` vs `zeros`
c1 = y > (th - hyst); % change #2 - precompute c1, c2
c2 = y > (th + hyst);
for k = 2:numel(y)
if yth(k - 1)
yth(k) = c1(k);
else
yth(k) = c2(k);
end
end
end
function yth = idea2()
c = [y > (th - hyst), y > (th + hyst)];
% Using Run-length encoding:
[v1,l1] = rlEncode(c(:,1));
[v2,l2] = rlEncode(c(:,2));
rle = cat(3,[v1 l1 cumsum(l1)],[v2 l2 cumsum(l2)]);
% rle(:,:,2) is similar to rle(:,:,1) with small changes:
% - col1 is circshifted by 1 and
% - col2 has the 1st and last elements switched
newenc = reshape([rle(1:2:end,3,2), rle(1:2:end,3,1)].', [], 1);
yth = rlDecode(rle(:,1,2), [newenc(1); diff(newenc(1:end-1))]);
end
function yth = idea3()
% Approximation of yth, should be almost indistinguishable with high resolution
yth = [0 0 repmat(repelem([1,0],50), 1, ceil(numel(y)/(2*50)) )].';
% The amount of zeros at the beginning as well as the value "50" are problem-specific
% and need to be computed using signal-specific considerations
yth = yth(1:1000);
end
function yth = idea4()
% Another approximation
yth = circshift(y > th, 1);
% The value "1" is problem-specific
end
end % q48637952
%% RLE (de)compression:
function [vals, lens] = rlEncode(vec)
J = find(diff([vec(1)-1; vec(:)]));
vals = vec(J);
lens = diff([J; numel(vec)+1]);
end
function vec = rlDecode(vals, lens)
l = cumsum([ 1; lens(:) ]);
z = zeros(1, l(end)-1);
z(l(1:end-1)) = 1;
vec = vals(cumsum(z));
end