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.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 % !!!! 0100 FCFACTOR=0.8; % put corner at Nyquist * 0.8 0101 fc=fc*FCFACTOR; 0102 % !!!! 0103 L = 2*N*pqmax + 1; 0104 h = p*firls( L-1, [0 2*fc 2*fc 1], [1 1 0 0]).*kaiser(L,bta)' ; 0105 % h = p*fir1( L-1, 2*fc, kaiser(L,bta)) ; 0106 else 0107 L = p; 0108 h = ones(1,p); 0109 end 0110 end 0111 0112 Lhalf = (L-1)/2; 0113 isvect = any(size(x)==1); 0114 if isvect 0115 Lx = length(x); 0116 else 0117 Lx = size(x, 1); 0118 end 0119 0120 % Need to delay output so that downsampling by q hits center tap of filter. 0121 nz = floor(q-mod(Lhalf,q)); 0122 z = zeros(1,nz); 0123 h = [z h(:).']; % ensure that h is a row vector. 0124 Lhalf = Lhalf + nz; 0125 0126 % Number of samples removed from beginning of output sequence 0127 % to compensate for delay of linear phase filter: 0128 delay = floor(ceil(Lhalf)/q); 0129 0130 % Need to zero-pad so output length is exactly ceil(Lx*p/q). 0131 nz1 = 0; 0132 while ceil( ((Lx-1)*p+length(h)+nz1 )/q ) - delay < ceil(Lx*p/q) 0133 nz1 = nz1+1; 0134 end 0135 h = [h zeros(1,nz1)]; 0136 0137 % ---- HERE'S THE CALL TO UPFIRDN ---------------------------- 0138 y = upfirdn(x,h,p,q); 0139 0140 % Get rid of trailing and leading data so input and output signals line up 0141 % temporally: 0142 Ly = ceil(Lx*p/q); % output length 0143 % Ly = floor((Lx-1)*p/q+1); <-- alternately, to prevent "running-off" the 0144 % data (extrapolation) 0145 if isvect 0146 y(1:delay) = []; 0147 y(Ly+1:end) = []; 0148 else 0149 y(1:delay,:) = []; 0150 y(Ly+1:end,:) = []; 0151 end 0152 0153 h([1:nz (end-nz1+1):end]) = []; % get rid of leading and trailing zeros 0154 % in case filter is output 0155 0156 0157 if numel(sz)>2 0158 y=reshape(y,[size(y,1),sz(2:end)]); 0159 end 0160