%% Moving Standard Deviation of Noisy Square Wave Signal
% Compute the moving standard deviation of a noisy square wave signal with 
% varying amplitude using the |dsp.MovingStandardDeviation| object.

%% Initialization
% Set up |movstdWindow| and |movstdExp| objects. 
% |movstdWindow| uses the sliding window method with a window length 
% of 800. |movstdExp| uses the exponential weighting method with a 
% forgetting factor of 0.999. Create a time scope for viewing the output.
FrameLength = 100;
Fs = 100;
movstdWindow = dsp.MovingStandardDeviation(800);
movstdExp = dsp.MovingStandardDeviation('Method','Exponential weighting',...
    'ForgettingFactor',0.999);
scope  = dsp.TimeScope('SampleRate',Fs,...
    'TimeSpanOverrunAction','Scroll',...
    'TimeSpan',1000,...
    'ShowGrid',true,...
    'BufferLength',1e7,...
    'YLimits',[0 3e-2]);
title = 'Actual (yellow) Sliding Window (blue) Exponentially Weighted (red) standard deviation';
scope.Title = title;

%% Compute the Standard Deviation
% Generate a noisy square wave signal. Vary the amplitude of the square wave 
% after a given number of frames. Apply the sliding window method and the 
% exponential weighting method on this signal. The actual standard deviation
% is |sqrt(np)|. This value is used while adding noise to the data. 
% Compare the actual standard deviation with the computed standard deviation 
% on the time scope.
count = 1;
noisepower = 1e-4 * [1 2 3 4];
index = 1;
for index = 1:length(noisepower)
    np = noisepower(index);
    yexp = sqrt(np)*ones(FrameLength,1);
    for i = 1:250
        x = sqrt(np) * randn(FrameLength,1);
        y1 = movstdWindow(x);
        y2 = movstdExp(x);
        scope([yexp,y1,y2]);
    end
end


%% 
% Copyright 2012 The MathWorks, Inc.