Home > NoiseTools > nt_resample.m

# nt_resample

## PURPOSE

RESAMPLE Change the sampling rate of a signal.

## SYNOPSIS

function [y, h] = resample( x, p, q, N, bta )

## DESCRIPTION

    RESAMPLE  Change the sampling rate of a signal.

    !!!!!
     Lower LP corner than Matlab's resample (Matlab's gives inadequate
     antialias filtering).
     Works with multidimensional arrays (Matlab's doesn't).
    !!!!!

       Y = RESAMPLE(X,P,Q) resamples the sequence in vector X at P/Q times
       the original sample rate using a polyphase implementation.  Y is P/Q
       times the length of X (or the ceiling of this if P/Q is not an integer).
       P and Q must be positive integers.

       RESAMPLE applies an anti-aliasing (lowpass) FIR filter to X during the
       resampling process, and compensates for the filter's delay.  The filter
       is designed using FIRLS.  RESAMPLE provides an easy-to-use alternative
       to UPFIRDN, relieving the user of the need to supply a filter or
       compensate for the signal delay introduced by filtering.

       In its filtering process, RESAMPLE assumes the samples at times before
       and after the given samples in X are equal to zero. Thus large
       deviations from zero at the end points of the sequence X can cause
       inaccuracies in Y at its end points.

       Y = RESAMPLE(X,P,Q,N) uses a weighted sum of 2*N*max(1,Q/P) samples of X
       to compute each sample of Y.  The length of the FIR filter RESAMPLE applies
       is proportional to N; by increasing N you will get better accuracy at the
       expense of a longer computation time.  If you don't specify N, RESAMPLE uses
       N = 10 by default.  If you let N = 0, RESAMPLE performs a nearest
       neighbor interpolation; that is, the output Y(n) is X(round((n-1)*Q/P)+1)
       ( Y(n) = 0 if round((n-1)*Q/P)+1 > length(X) ).

       Y = RESAMPLE(X,P,Q,N,BTA) uses BTA as the BETA design parameter for the
       Kaiser window used to design the filter.  RESAMPLE uses BTA = 5 if
       you don't specify a value.

       Y = RESAMPLE(X,P,Q,B) uses B to filter X (after upsampling) if B is a
       vector of filter coefficients.  RESAMPLE assumes B has odd length and
       linear phase when compensating for the filter's delay; for even length
       filters, the delay is overcompensated by 1/2 sample.  For non-linear
       phase filters consider using UPFIRDN.

       [Y,B] = RESAMPLE(X,P,Q,...) returns in B the coefficients of the filter
       applied to X during the resampling process (after upsampling).

       If X is a matrix, RESAMPLE resamples the columns of X.

       % Example:
       %   Resample a sinusoid at 3/2 the original rate.

       tx = 3:3:300;           % Time vector for original signal
       x = sin(2*pi*tx/300);   % Define a sinusoid
       ty = 2:2:300;           % Time vector for resampled signal
       y = resample(x,3,2);    % Change sampling rate
       plot(tx,x,'+-',ty,y,'o:')
       legend('original','resampled'); xlabel('Time')

       See also UPFIRDN, INTERP, DECIMATE, FIRLS, KAISER, INTFILT.

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 function  [y, h] = resample( x, p, q, N, bta )
    0002 %RESAMPLE  Change the sampling rate of a signal.
    0003 %
    0004 %!!!!!
    0005 % Lower LP corner than Matlab's resample (Matlab's gives inadequate
    0006 % antialias filtering).
    0007 % Works with multidimensional arrays (Matlab's doesn't).
    0008 %!!!!!
    0009 %
    0010 %   Y = RESAMPLE(X,P,Q) resamples the sequence in vector X at P/Q times
    0011 %   the original sample rate using a polyphase implementation.  Y is P/Q
    0012 %   times the length of X (or the ceiling of this if P/Q is not an integer).
    0013 %   P and Q must be positive integers.
    0014 %
    0015 %   RESAMPLE applies an anti-aliasing (lowpass) FIR filter to X during the
    0016 %   resampling process, and compensates for the filter's delay.  The filter
    0017 %   is designed using FIRLS.  RESAMPLE provides an easy-to-use alternative
    0018 %   to UPFIRDN, relieving the user of the need to supply a filter or
    0019 %   compensate for the signal delay introduced by filtering.
    0020 %
    0021 %   In its filtering process, RESAMPLE assumes the samples at times before
    0022 %   and after the given samples in X are equal to zero. Thus large
    0023 %   deviations from zero at the end points of the sequence X can cause
    0024 %   inaccuracies in Y at its end points.
    0025 %
    0026 %   Y = RESAMPLE(X,P,Q,N) uses a weighted sum of 2*N*max(1,Q/P) samples of X
    0027 %   to compute each sample of Y.  The length of the FIR filter RESAMPLE applies
    0028 %   is proportional to N; by increasing N you will get better accuracy at the
    0029 %   expense of a longer computation time.  If you don't specify N, RESAMPLE uses
    0030 %   N = 10 by default.  If you let N = 0, RESAMPLE performs a nearest
    0031 %   neighbor interpolation; that is, the output Y(n) is X(round((n-1)*Q/P)+1)
    0032 %   ( Y(n) = 0 if round((n-1)*Q/P)+1 > length(X) ).
    0033 %
    0034 %   Y = RESAMPLE(X,P,Q,N,BTA) uses BTA as the BETA design parameter for the
    0035 %   Kaiser window used to design the filter.  RESAMPLE uses BTA = 5 if
    0036 %   you don't specify a value.
    0037 %
    0038 %   Y = RESAMPLE(X,P,Q,B) uses B to filter X (after upsampling) if B is a
    0039 %   vector of filter coefficients.  RESAMPLE assumes B has odd length and
    0040 %   linear phase when compensating for the filter's delay; for even length
    0041 %   filters, the delay is overcompensated by 1/2 sample.  For non-linear
    0042 %   phase filters consider using UPFIRDN.
    0043 %
    0044 %   [Y,B] = RESAMPLE(X,P,Q,...) returns in B the coefficients of the filter
    0045 %   applied to X during the resampling process (after upsampling).
    0046 %
    0047 %   If X is a matrix, RESAMPLE resamples the columns of X.
    0048 %
    0049 %   % Example:
    0050 %   %   Resample a sinusoid at 3/2 the original rate.
    0051 %
    0052 %   tx = 3:3:300;           % Time vector for original signal
    0053 %   x = sin(2*pi*tx/300);   % Define a sinusoid
    0054 %   ty = 2:2:300;           % Time vector for resampled signal
    0055 %   y = resample(x,3,2);    % Change sampling rate
    0056 %   plot(tx,x,'+-',ty,y,'o:')
    0057 %   legend('original','resampled'); xlabel('Time')
    0058 %
    0059 %   See also UPFIRDN, INTERP, DECIMATE, FIRLS, KAISER, INTFILT.
    0060
    0061 %   NOTE-1: digital anti-alias filter is desiged via windowing
    0062
    0063 %   Author(s): James McClellan, 6-11-93
    0064 %              Modified to use upfirdn, T. Krauss, 2-27-96
    0065 %   Copyright 1988-2011 The MathWorks, Inc.
    0066 %
    0067
    0068 if nargin < 5,  bta = 5;  end   %--- design parameter for Kaiser window LPF
    0069 if nargin < 4,   N = 10;   end
    0070 if abs(round(p))~=p || p==0
    0071   error(message('signal:resample:MustBePosInteger', 'P'));
    0072 end
    0073 if abs(round(q))~=q || q==0
    0074   error(message('signal:resample:MustBePosInteger', 'Q'));
    0075 end
    0076
    0077 sz=size(x);
    0078 if numel(sz)>2
    0079     x=reshape(x,sz(1),prod(sz(2:end)));
    0080 end
    0081
    0082 [p,q] = rat( p/q, 1e-12 );  %--- reduce to lowest terms
    0083    % (usually exact, sometimes not; loses at most 1 second every 10^12 seconds)
    0084 if (p==1) && (q==1)
    0085     y = x;
    0086     h = 1;
    0087     if numel(sz)>2
    0088         y=reshape(y,[size(y,1),sz(2:end)]);
    0089     end
    0090     return
    0091 end
    0092 pqmax = max(p,q);
    0093 if length(N)>1      % use input filter
    0094    L = length(N);
    0095    h = N;
    0096 else                % design filter
    0097    if( N>0 )
    0098       fc = 1/2/pqmax;
    0099       FCFACTOR=0.8; % put corner at Nyquist * 0.8
    0100       fc=fc*FCFACTOR;
    0101       L = 2*N*pqmax + 1;
    0102       h = p*firls( L-1, [0 2*fc 2*fc 1], [1 1 0 0]).*kaiser(L,bta)' ;
    0103       % h = p*fir1( L-1, 2*fc, kaiser(L,bta)) ;
    0104    else
    0105       L = p;
    0106       h = ones(1,p);
    0107    end
    0108 end
    0109
    0110 Lhalf = (L-1)/2;
    0111 isvect = any(size(x)==1);
    0112 if isvect
    0113     Lx = length(x);
    0114 else
    0115     Lx = size(x, 1);
    0116 end
    0117
    0118 % Need to delay output so that downsampling by q hits center tap of filter.
    0119 nz = floor(q-mod(Lhalf,q));
    0120 z = zeros(1,nz);
    0121 h = [z h(:).'];  % ensure that h is a row vector.
    0122 Lhalf = Lhalf + nz;
    0123
    0124 % Number of samples removed from beginning of output sequence
    0125 % to compensate for delay of linear phase filter:
    0126 delay = floor(ceil(Lhalf)/q);
    0127
    0128 % Need to zero-pad so output length is exactly ceil(Lx*p/q).
    0129 nz1 = 0;
    0130 while ceil( ((Lx-1)*p+length(h)+nz1 )/q ) - delay < ceil(Lx*p/q)
    0131     nz1 = nz1+1;
    0132 end
    0133 h = [h zeros(1,nz1)];
    0134
    0135 % ----  HERE'S THE CALL TO UPFIRDN  ----------------------------
    0136 y = upfirdn(x,h,p,q);
    0137
    0138 % Get rid of trailing and leading data so input and output signals line up
    0139 % temporally:
    0140 Ly = ceil(Lx*p/q);  % output length
    0141 % Ly = floor((Lx-1)*p/q+1);  <-- alternately, to prevent "running-off" the
    0142 %                                data (extrapolation)
    0143 if isvect
    0144     y(1:delay) = [];
    0145     y(Ly+1:end) = [];
    0146 else
    0147     y(1:delay,:) = [];
    0148     y(Ly+1:end,:) = [];
    0149 end
    0150
    0151 h([1:nz (end-nz1+1):end]) = [];  % get rid of leading and trailing zeros
    0152                                  % in case filter is output
    0153
    0154
    0155 if numel(sz)>2
    0156     y=reshape(y,[size(y,1),sz(2:end)]);
    0157 end
    0158