middlepad.m

linear time－frequency toolbox

function f=middlepad(f,L,centering)
%MIDDLEPAD  Symmetrically zero-extends or cuts a function.
%   Usage:  h=middlepad(f,L);
%
%   MIDDLEPAD(f,L) cuts or zero-extends f to length L by inserting
%   zeros in the middle of the vector, or by cutting in the middle
%   of the vector.
%
%   If f is whole-point even, MIDDLEPAD(f,L) will also be whole-point even.
%
%   If f has even length, then f will not be purely zero-extended, but
%   the last element will be repeated once and multiplied by 1/2!
%   That is, the support of f will increase by one!
%
%   MIDDLEPAD can be used to Fourier interpolate a vector. The following
%   code will return an interpolation of f to length L:
% 
%       f=idft(middlepad(dft(f),L));
% 
%   SEE ALSO:  ISEVEN, FIREXTEND

error(nargchk(2,3,nargin));

if nargin==2
  centering=0;
end;

if  (prod(size(L))~=1 || ~isnumeric(L))
  error([callfun,': L must be a scalar']);
end;

if L<1
  error('L must be larger than 0.');
end;

Lorig=length(f);

if Lorig==L
  % Nothing to do.
  return;
end;

wasrow=0;
if size(f,2)>1
  if size(f,1)>1
    error('f must be a vector');
  else
    % f was a row vector.
    wasrow=1;
    f=f(:);
  end;
end;

switch centering
  case 0

    % ---------------   WPE case --------------------------------------
    
    if Lorig==1
      % Rather trivial case
      f=[f(1);zeros(L-1,1)];
      
    else
      if Lorig>L
	% Cut
	
	if mod(L,2)==0
	  
	  % L even. Use average of endpoints.
	  f=[f(1:L/2);(f(L/2+1)+f(Lorig-L/2+1))/2;f(Lorig-L/2+2:Lorig)];
	  
	else
	  
	  % No problem, just cut.
	  f=[f(1:(L+1)/2);f(Lorig-(L-1)/2+1:Lorig)];
	  
	end;     
	
      else
	
	d=L-Lorig;
	
	% Extend
	if mod(Lorig,2)==0
	  
	  % Lorig even. We must split a value.
	  
	  f=[f(1:Lorig/2);...
	     f(Lorig/2+1)/2;...
	     zeros(d-1,1);...
	     f(Lorig/2+1)/2;...
	     f(Lorig/2+2:Lorig)];
	  
	else
	  % Lorig is odd, we can just insert zeros.
	  f=[f(1:(Lorig+1)/2);zeros(d,1);f((Lorig+3)/2:Lorig)];
	  
	end;
	
      end;
    end;
    
  case .5

    % ------------------ HPE case ------------------------------------
    
    if Lorig==1
      
    else
      if Lorig>L

	d=Lorig-L;
	% Cut
	
	if mod(L,2)==0
	  % L even
	  
	  % No problem, just cut.
	  f=[f(1:L/2);...
	     f(Lorig-L/2+1:Lorig);];

	else

	  % Average of endpoints.
	  f=[f(1:(L-1)/2);(f((L+1)/2)+f(Lorig-(L-1)/2))/2;...
	     f(Lorig-(L-1)/2+1:Lorig);];
	    
	end;
		  
      else
	
	d=L-Lorig;
	
	% Extend
	if mod(Lorig,2)==0
	  
	  % Lorig even. We can just insert zeros in the middle.
	  
	  f=[f(1:Lorig/2);...
	     zeros(d,1);...
	     f(Lorig/2+1:Lorig)];
	  
	else
	  % Lorig odd. We need to split a value in two
	  f=[f(1:(Lorig-1)/2);...
	     f((Lorig+1)/2)/2;...
	     zeros(d-1,1);...
	     f((Lorig+1)/2)/2;...
	     f((Lorig-1)/2+2:Lorig)];
	  
	end;
	
      end;

    end;

  otherwise
    error('Unsupported centering.');
end;

if wasrow
  f=f.';
end;