Berlekamp-Massey in Matlab

So I'm trying to implement the inversionless Berlekamp-Massey algorithm in
Matlab, and I'm having a devil of a time.  I've gotten it to the point that
it'll calculate correct polynomials for a simple (7,3) code (from a chapter
by Bernard Sklar on the subject) but it falls apart when I give it a larger
(204,188) code (from ITU-T J.83), and I have no idea why. 

Could it have something to do with the fact the J.83 code uses roots 0-15
instead of 1-16, does BM care?

At any rate, I figure I have nothing to lose at this point by just posting
the code I have thus far and seeing if anyone can immediately spot what I'm
doing wrong.  


function [loc_poly,eval_poly] = rs_berlekamp(m,pp,syndrome)
% function [loc_poly,eval_poly] = rs_berlekamp(syndrome)
%   Finds the error locator and error evaluator polynomials given the
%   syndrome from a received codeword
%
%   m        - Number of bits in each symbol, sets the field to GF(2^m)
%   pp       - Primitive polynomial for the field
%   syndrome - Syndrome computed from received codeword 

% Calculate size of syndrome and from that our error correction
capability.
s_len = length(syndrome);
T     = s_len/2;

% Initialize values
lambda = gf(horzcat([1], zeros(1,T)), m, pp);
beta   = lambda;

kappa  = 0;
gamma  = gf(1,m,pp);  % Must be a field element b/c we'll multiply by it

% Loop for 2T iterations
for r = 1:(2*T)
    % Step 1, convolve syndrome and previous lambda to get delta
    delta = gf(0,m,pp);
    for ii = 0:min(T,r-1)
        delta = delta + (syndrome(r-ii)*lambda(ii+1));
    end
    
    % Step 2, update lambda value
    lambda(1) = gamma*lambda(1);
    for ii = 1:T
        lambda(ii+1) = gamma*lambda(ii+1)-delta*beta(ii);
    end
    
    % Step 3, update variables
    if ((delta ~= gf(0,m,pp)) && (2*kappa <= r-1))
        beta  = lambda;
        gamma = delta;
        kappa = r-kappa;
    else
        % Multiply beta by X (equivalent to a shift)
        beta(2:T+1) = beta(1:T);
        beta(1) = 0;
    end
end

loc_poly = lambda;

% Step 4, convolve syndrome and loc_poly to get error evaluator
polynomial
eval_poly = gf(zeros(1,T),m,pp);

for ii=0:T-1
    eval_poly(ii+1) = gf(0,m,pp);
    for jj=0:ii
        eval_poly(ii+1) = eval_poly(ii+1) +
syndrome(ii-jj+1)*loc_poly(ii+1);
    end
end


Some test vectors:
For the (7,3) code a syndrome of [0 5 7 3] should yield a locator
polynomial with roots at a^3 and a^4 (a==2, default polynomial in matlab
for GF(8)).  The error valuator run through forney's algorithms should give
a^5 and a^0 for error values (7 and 1)

For the J.83 code, I used this test vector:

tv_1 = gf([76 103 77 193 92 32 157 91 93 17 221 117 20 231 72 157 169 51
245 170 138 222 142 5 123 206 188 146 2 183 115 168 192 205 7 199 145 19 64
34 144 138 18 252 64 80 77 11 135 65 104 242 234 31 151 92 183 134 67 126
218 185 51 40 94 220 175 162 36 20 223 107 124 117 131 69 59 229 232 154 93
153 170 228 22 137 109 157 143 58 27 3 15 82 199 86 158 253 165 138 59 189
227 219 152 167 233 110 74 161 75 159 12 254 53 155 89 183 7 17 236 22 85
134 63 138 199 133 238 38 106 71 153 9 16 82 25 43 95 10 181 164 44 16 104
118 52 222 152 6 229 116 15 27 255 221 157 7 82 118 25 146 83 115 147 19 15
77 133 143 62 233 211 113 250 147 60 207 115 64 244 36 131 248 165 157 120
147],8);

Which, during encoding, has the following appended:
[63 70 180 135 188 126 6 130 178 247 43 46 143 228 136 44]

I found the syndrome of the non-corrupted codeword and it did indeed come
out to zero like it should.  So I took the second element (103) and
simulated a symbol error by setting it to zero.  I should get a locator
polynomial with a root at a^1 correct?  

I should note that I had to flip the syndrome for the (7,3) code from [3 7
5 0] before it would work, I don't know if that's diagnostic or not.

On 2009-03-19, gct <smcallis@gmail.com> wrote:
> by Bernard Sklar on the subject) but it falls apart when I give it a larger
> (204,188) code (from ITU-T J.83)

You do realize that's a "shortened" (zero-padded) 255,239 code, right?
Also that only applies to Annex A and C which are not used in North
America.  We use J.83 Annex B which is much uglier overall and does the
RS step over 7 bit symbols.

-- 
Ben Jackson AD7GD
<ben@ben.com>
http://www.ben.com/

>On 2009-03-19, gct <smcallis@gmail.com> wrote:
>> by Bernard Sklar on the subject) but it falls apart when I give it a
larger
>> (204,188) code (from ITU-T J.83)
>
>You do realize that's a "shortened" (zero-padded) 255,239 code, right?
>Also that only applies to Annex A and C which are not used in North
>America.  We use J.83 Annex B which is much uglier overall and does the
>RS step over 7 bit symbols.
>
>-- 
>Ben Jackson AD7GD
><ben@ben.com>
>http://www.ben.com/
>

Sure enough, I'm not actually building a J.83 modem I'm just using it as a
reasonable "real world" example to test out this implementation.  I
understand it's shortened.  It was my understand that preprending zeros
like the spec says doesn't really affect anything.  Would it change how you
decode it?

On 2009-03-19, gct <smcallis@gmail.com> wrote:
>>> (204,188) code (from ITU-T J.83)
>>
>>You do realize that's a "shortened" (zero-padded) 255,239 code, right?
>
> It was my understand that preprending zeros
> like the spec says doesn't really affect anything.  Would it change how you
> decode it?

I believe you append the zeroes.  You have to push them through the input
to get all the state out on the encode side.  On the receive side you can
discard them.  Knowing they're zero gives you information when looking for
errors.

I'm not familiar with the algorithm you're trying to implement in Matlab;
I was just wondering if you were putting in the right parameters.

-- 
Ben Jackson AD7GD
<ben@ben.com>
http://www.ben.com/