off by two minimum phase signal from 'rceps' for odd length input

[Sergei Steshenko]

Hello,

here is a quick example:

"
octave:43> [y, xm] = rceps([1 2 3 4]);
octave:44> xm
xm =

   3.98049   3.01913   2.01951   0.98087

octave:45> [y, xm] = rceps([1 2 3 4 5]);
octave:46> xm
xm =

   4.9244   3.7654   1.5151

octave:47>     

".

In the second case 'xm' has just 3 elements, though it is supposed to have 5.

Source code of 'rceps' explicitly separately deals with odd and even length 
inputs; I suspect that

n/2-1

in "zeros(1,n/2-1)", "zeros(n/2-1,columns(y))" pieces in places related to odd 
length input should be replaced with

n/2+1

, but I am not sure. I.e. I am not sure what to do with the components around 
Nyquist frequency for odd length input.


So, is this off by two signal a bug to bug compatibility with Matlab or just a 
bug ?

Thanks,
  Sergei.

P.S. "Fixed" by me as described above code without the comments and test cases:

function [y, ym] = rceps(x)
  fprintf(stderr, "my_rceps: numel(x)=%u\n", numel(x));

  if (nargin != 1)
    print_usage;
  end
  y = real(ifft(log(abs(fft(x)))));
  if nargout == 2
    n=length(x);
    if rows(x)==1
      if rem(n,2)==1
        ym = [y(1), 2*y(2:n/2), zeros(1,n/2+1)]; # changed line
      else
        ym = [y(1), 2*y(2:n/2), y(n/2+1), zeros(1,n/2-1)];
      endif
    else
      if rem(n,2)==1
        ym = [y(1,:); 2*y(2:n/2,:); zeros(n/2+1,columns(y))]; # changed line
      else
        ym = [y(1,:); 2*y(2:n/2,:); y(n/2+1,:); zeros(n/2-1,columns(y))];
      endif
    endif
    fprintf(stderr, "my_rceps: numel(ym)=%u\n", numel(ym));
    ym = real(ifft(exp(fft(ym))));
    fprintf(stderr, "my_rceps: numel(ym)=%u\n", numel(ym));
  endif
endfunction


Original code:

     33 function [y, ym] = rceps(x)
     34   if (nargin != 1)
     35     print_usage;
     36   end
     37   y = real(ifft(log(abs(fft(x)))));
     38   if nargout == 2
     39     n=length(x);
     40     if rows(x)==1
     41       if rem(n,2)==1
     42         ym = [y(1), 2*y(2:n/2), zeros(1,n/2-1)];
     43       else
     44         ym = [y(1), 2*y(2:n/2), y(n/2+1), zeros(1,n/2-1)];
     45       endif
     46     else
     47       if rem(n,2)==1
     48         ym = [y(1,:); 2*y(2:n/2,:); zeros(n/2-1,columns(y))];
     49       else
     50         ym = [y(1,:); 2*y(2:n/2,:); y(n/2+1,:); 
zeros(n/2-1,columns(y))];
     51       endif
     52     endif
     53     ym = real(ifft(exp(fft(ym))));
     54   endif
     55 endfunction