How does this Matlab code work?

Question

I came across this uncommented Matlab code except for the header:

*************************************
function y=movingvar(X,N)
% y=movingvar(X,N)
% Calculates N-point moving variance of Vector X
% Highly recommend that N be odd (no error checking)
% Note: first and last N/2 points will be unreliable.
% Output will be a column vector.
% Authors: Scott Seidman ([email protected]) 1/23/99

X=X(:);
XSQR=X.*X;
convsig=ones(1,N);
y=(conv(convsig,XSQR)-(conv(convsig,X).^2)/N)/(N-1);

y=y(ceil(N/2):length(X)+floor(N/2));
****************************************************

(from http://www.mathworks.com/matlabcentral/newsreader/view_thread/39521)

According to the comments here, it can also be used as a realtime RMS calculator: https://electronics.stackexchange.com/a/135276/53375

I can see that being useful, but I don't understand Matlab. Could someone please explain what's going on?

I'm academically familiar with convolution, and I do know C.

Answer 1

Starting with the formula for sample variance $$ s^{2}=\frac{ \sum_{i=1}^{N}\left ( X_{i}-\bar{X} \right )^2}{N-1} $$

you multiply out to get $$ s^{2}= \frac{\sum_{i=1}^{N}\left({X_i}^2-2X_i \bar{X}+\bar{X}^2\right )}{N-1} $$

Distribute the sums. Using the fact that \$ \bar{X}\$ is a constant, then $$ \sum 2X_i \bar{X} = 2\bar{X}\sum{X_i}, $$ and $$\sum{X_i}=N\bar{X},$$

So, $$ 2\bar{X}\sum{X_i}=2N\bar{X}^2.$$

Then, using $$ \sum_{i=1}^{N}\bar{X}^2=N\bar{X}^2, $$ and $$ N\bar{X}^2 = \frac{ N\left( \sum{X} \right )^2}{N^2} $$

the overall equation for varance reduces to a convenient computational formula for variance: $$ \frac{ \sum_{i-1}^{N}\left ( X_{i}^2 \right )-\frac{\left (\sum_{i=1}^{N}X_i \right )^2}{N}}{N-1} $$

XSQR=X.*X; \ simply element by element squares the column array, X

Without entirely explaining convolution, convolution with a string of ones of length N (i.e., convsig) is simply a sum of the previous N points at each time step.

The last line in the code just gets rid of the trash at either end of the array associated with the beginning and end of convolution.

Lastly, you need to understand that most operations on arrays in Matlab operate on the whole array at the same time -- most of the operations in the program you show would absolutely need to loop over every element in the array in a language like C, which likely very confusing for folks not used to working in Matlab.

Bottom line, its not much different than actually just implementing the sum, but convolution is something directly written into matlab, and there are generally libraries that directly handle convolution on DSP's, so its a convenient way to implement a moving variance filter without needing to do the grunt work.