function y2 = decim_eval(N, M)

x = [1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0];
x = randn(1, 1000);
x = 1:20;
w = FIRCalcLowpass(0.5, N);
w = 1:N;

y = filter(w, 1, x)';
y_dec = y(1:M:lge(y))
y_dec2 = decim(w, M, N, x)'

y2 = y_dec2 - y_dec;

function [y] = decim(w, M, N, x)
Nout = lge(x(1:M:lge(x)));

for k=0:M-1
    nz = fix(k/M) + (k > 0);
    d = (M-k) * (k > 0);
    xp = [zeros(1, nz) x(d+1:M:lge(x))]
%    xp = [zeros(1, k) x];
%    xp = xp(1:M:lge(x))
    wp = w(k+1:M:N);
    yp = filter(wp, 1, xp);
    summer(k+1, :) = yp(1:Nout);
end;
y = sum(summer, 1);

return;