Search code examples
matlabinterpolationenvelope

Encase signal by getting the upper and lower envelope in MATLAB

I have a signal that bounces around some values. I want to encase the signal from above and below as seen in the following figure: example

If the signal only grows or decreases, using a running maxima script works great:

function [out] = rMax(X)
    Y=X;
    maximum=X(1);
    for k=1:length(X)
        if X(k)<=maximum
        Y(k)=maximum;
        else
            maximum=X(k);
        end
    end
    out=Y;
end

However, when the signal alternates, I cannot use that method anymore. Is there a way to achieve this in MATLAB or Mathematica?

Solution

  • You can use movmax and movmin in conjunction with appropriate window sizes n1 > n2 depending on your needs. These two functions were introduced in R2016a. Have a look at the bottom of my answer for a self-written replacement for them that works prior to R2016a.

    To get the upper bounds, you can use:

    xMax = movmax(x,[n1,n1]);
xMax = movmin(xMax,[n2,n2]);

    For the lower bounds, just switch movmin and movmax:

    xMin = movmin(x,[n1,n1]);
xMin = movmax(xMin,[n2,n2]);

    The result can look like this:

    result

    If you choose n2=n1, the bounds are very tight when there is a peak in the data x. In case you choose a bigger difference by just making n2 smaller than n1, you will get a longer straight line at the peak. A side effect is that the bounds start to be distanced when the signal jumps from 1 to -1 or vice versa. If n1 is chosen too high, the small part of the signal around 100 won't be extracted.

    Here is the complete code to generate the above figure with some sample data matching the plot in your question. Just play around with the values n1 and n2 to see their effect.

    n1 = 20;        % for first window
n2 = 18;        % for second window

% generate sample data
t = 20:0.1:220;
x = -ones(size(t));
x(t>60&t<100) = 1;
x(t>105&t<135) = 1;
x(t>145&t<155) = 1;
x = x + 0.4*randn(size(x));

% get upper bounds
xMax = movmax(x,[n1,n1]);
xMax = movmin(xMax,[n2,n2]);

% get lower bounds
xMin = movmin(x,[n1,n1]);
xMin = movmax(xMin,[n2,n2]);

% draw figure for illustration
figure; hold on;
plot(t,x);
xlim([20,220]);
ylim([-3,3]);
plot(t,xMax,'r','LineWidth',1.1);
plot(t,xMin,'Color',[0,0.5,0],'LineWidth',1.1);

    Prior to R2016a

    To have the basic functionality of movmin and movmax in MATLAB versions prior to R2016a, we can implement our own functions. Therefore, we need to apply a moving-minimum (or moving-maximum respectively) which can be achieved quite easily. In order to maintain compatibility with the functions in R2016a, we are going to implement the case where the second argument is a scalar and where it is a vector of two elements. This covers the following syntax similar as seen here with the restriction that x needs to be a vector.

    y = movingmax(x,k)
y = movingmax(x,[kb kf])

    Here is the code for movingmax that replaces movmax:

    function y = movingmax(x,n)
if length(n) > 1
    a=n(1); b=n(2);
else
    b=floor((n-1)/2); a=b+mod(n-1,2);
end
s = size(x);
xp = [ones(a,1)*x(1);x(:);ones(b,1)*x(end)];
y = zeros(size(x));
for k = 1:length(x)
    y(k) = max(xp(k:k+a+b));
end
y = reshape(y,s);

    Here is the code for movingmin that replaces movmin:

    function y = movingmin(x,n)
if length(n) > 1
    a=n(1); b=n(2);
else
    b=floor((n-1)/2); a=b+mod(n-1,2);
end
s = size(x);
xp = [ones(a,1)*x(1);x(:);ones(b,1)*x(end)];
y = zeros(size(x));
for k = 1:length(x)
    y(k) = min(xp(k:k+a+b));
end
y = reshape(y,s);