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

Generated on Mon 10-Nov-2014 14:40:42 by m2html © 2005