Documentation of idresamp

Global Index (short | long) | Local contents | Local Index (short | long)

Function Synopsis

[zd,decr] = idresamp(zdat,r,nfilt,tol)

Help text

IDRESAMP Resamples data by decimation and interpolation.
   ZD = IDRESAMP(Z,R)

   Z  : The output-input data matrix Z = [Y U]. Multi-variate
        data are allowed
   R  : The resampling factor. The new data record ZD will correspond
        to a sampling interval of R times that of the original data.
        R > 1 thus means decimation and R < 1 means interpolation.
        Any positive number for R is allowed, but it will be replaced
        by a rational approximation.
   ZD : The resampled data. The columns of ZD correspond to those of Z.

   [ZD, ACT_R] = IDRESAMP(Z,R,ORDER,TOL) gives access to the following
   options: ORDER determines the filter orders used at decimation
   and interpolation (Default 8). TOL gives the tolerance of the
   rational approximation (Default 0.1). ACT_R is the actually used
   resampling factor.
   See also IDFILT.

Cross-Reference Information

Listing of function idresamp

function [zd,decr] = idresamp(zdat,r,nfilt,tol)

%   L. Ljung 3-3-95.
%   The code is adopted from the routines DECIMATE and INTERP
%   in the Signal Processing Toolbox, written by L. Shure.
%   Copyright (c) 1986-98 by The MathWorks, Inc.
%   $Revision: 3.3 $  $Date: 1997/12/02 03:44:29 $

if nargin < 2
   disp('Usage: ZD = IDRESAMP(Z,R).')
   disp('       [ZD,ACT_R] = IDRESAMP(Z,R,FILTER_ORDER,TOLERANCE)')
   return
end

if nargin<4
   tol=[];
end
if nargin<3
   nfilt=[];
end

if isempty(tol),tol=0.1;end
if isempty(nfilt),nfilt=8;end
if r<=0, error('The resampling factor R must be positive.'),end
[ndat,nyu]=size(zdat);

if ndat<nyu
   error('Data should be organized columnwise.')
end
if ndat<10
 error('The data record length must be at least 10.')
end
[numr,denr] = rat(r,tol);
decr = numr/denr;
if numr == 1 & denr == 1
        zd = zdat;
        return
end
if denr>1

   l = nfilt/2;
   alpha = .5;
   % calculate AP and AM matrices for inversion
   s1 = toeplitz(0:l-1) + eps;
   s2 = hankel(2*l-1:-1:l);
   s2p = hankel([1:l-1 0]);
   s2 = s2 + eps + s2p(l:-1:1,l:-1:1);
   s1 = sin(alpha*pi*s1)./(alpha*pi*s1);
   s2 = sin(alpha*pi*s2)./(alpha*pi*s2);
   ap = s1 + s2;
   am = s1 - s2;
   ap = inv(ap);
   am = inv(am);

   % now calculate D based on INV(AM) and INV(AP)
   d = zeros(2*l,l);
   d(1:2:2*l-1,:) = ap + am;
   d(2:2:2*l,:) = ap - am;
   % set up arrays to calculate interpolating filter B
   x = (0:denr-1)/denr;
   y = zeros(2*l,1);
   y(1:2:2*l-1) = (l:-1:1);
   y(2:2:2*l) = (l-1:-1:0);
   X = ones(2*l,1);
   X(1:2:2*l-1) = -ones(l,1);
   XX = eps + y*ones(1,denr) + X*x;
   y = X + y + eps;
   h = .5*d'*(sin(pi*alpha*XX)./(alpha*pi*XX));
   b = zeros(2*l*denr+1,1);
   b(1:l*denr) = h';
   b(l*denr+1) = .5*d(:,l)'*(sin(pi*alpha*y)./(pi*alpha*y));
   b(l*denr+2:2*l*denr+1) = b(l*denr:-1:1);
   for kc = 1:nyu
       idata = zdat(:,kc);
       % use the filter B to perform the interpolation
       odata = zeros(denr*ndat,1);
       odata(1:denr:ndat*denr) = idata;
       % Filter a fabricated section of data first
       od = zeros(2*l*denr,1);
       od(1:denr:(2*l*denr)) = 2*idata(1)-idata((2*l+1):-1:2);
       [od,zi] = filter(b,1,od);
       [odata,zf] = filter(b,1,odata,zi);
       odata(1:(ndat-l)*denr) = odata(l*denr+1:ndat*denr);

       od = zeros(2*l*denr,1);
       od(1:denr:(2*l)*denr) = ...
                        [2*idata(ndat)-(idata((ndat-1):-1:(ndat-2*l)))];
       od = filter(b,1,od,zf);
       odata(ndat*denr-l*denr+1:ndat*denr) = od(1:l*denr);
       zd(:,kc) = odata;
    end
    zdat = zd;
    [ndat,dum]=size(zdat);
end
if numr>1
   nout = ceil(ndat/numr);
   odata = idfilt(zdat,nfilt,0.8/numr);
   nbeg = numr - (numr*nout - ndat);
   zd = odata(nbeg:numr:ndat,:);
end