
# my_resample.m -- Resample a time-series by a rational factor

# This file may be used under the terms defined by the GNU 
# General Public License [write to the Free Software Foundation, 
# 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA]
#
# WITHOUT ANY WARRANTY; without even the implied warranty 
# of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
#
# Copyright (C) 2008 Eric Chassande-Mottin, CNRS (France)

#

function  [y, h] = my_resample( x, p, q, h )

nargchk(3,4,nargin);

if any([p q]<=0)|any([p q]~=floor([p q])),
  error("resample.m: p and q must be positive integers");
endif;

## simplify decimation and interpolation factors

great_common_divisor=gcd(p,q);
if (great_common_divisor>1)
  p=p/great_common_divisor;
  q=q/great_common_divisor;
endif;

## filter design if required

if (nargin < 4)

 rejection=-3;
 stopband_cutoff_f = 1.0/(2.0 * max(p,q));
 roll_off_width = stopband_cutoff_f / 10.0;
 
 ## determine filter length
 ## use empirical formula from Oppenheim and Schafer, Discrete-Time Signal Proc. Prentice Hall. Chap 7, Eq. (7.63) p 476

 rejection_dB=-20.0*rejection;
 L=ceil((rejection_dB-8.0) / (28.714 * roll_off_width));
    
 ## ideal sinc filter

 t=(-L:L)';
 ideal_filter=2*p*stopband_cutoff_f*sinc(2*stopband_cutoff_f*t);  

 ## determine parameter of Kaiser window
 ## use empirical formula from Oppenheim and Schafer, Discrete-Time Signal Proc. Prentice Hall. Chap 7, Eq. (7.62) p 474

 if ((rejection_dB>=21)&(rejection_dB<=50))
   beta = 0.5842 * (rejection_dB-21.0)^0.4 + 0.07886 * (rejection_dB-21.0);
 elseif (rejection_dB>50)
   beta = 0.1102 * (rejection_dB-8.7);
 else
   beta = 0.0;
 endif;

 ## apodize ideal filter

 h=kaiser(2*L+1,beta).*ideal_filter;
 
endif;

Lx = length(x);
Lh = length(h);
Ly = ceil(Lx*p/q);

## pre and postpad filter response

nz_pre = floor(q-mod(L,q));
hpad = prepad(h,Lh+nz_pre);

offset = floor((L+nz_pre)/q);
nz_post = 0;
while ceil( ( (Lx-1)*p + nz_pre + Lh + nz_post )/q ) - offset < Ly
    nz_post++;
endwhile;
hpad = postpad(hpad,Lh + nz_pre + nz_post);

## filtering

[row col]=size(x);

if (row==1) || (col==1), # input is a vector
  xfilt = upfirdn(x,hpad,p,q);
  y = xfilt(offset+1:offset+Ly);
else # input is a matrix
  y = zeros(Ly,col);
  for n = 1:col,
    xfilt=upfirdn(x(:,n),hpad,p,q);
    y(:,n)=xfilt(offset+1:offset+Ly);
  endfor;
endif

