www.gusucode.com > signal 工具箱matlab源码程序 > signal/interp.m
function [odata,b] = interp(idata,r,l,cutoff) %INTERP Resample data at a higher rate using lowpass interpolation. % Y = INTERP(X,R) resamples the sequence in vector X at R times % the original sample rate. The resulting resampled vector Y is % R times longer, LENGTH(Y) = R*LENGTH(X). % % A symmetric filter, B, allows the original data to pass through % unchanged and interpolates between so that the mean square error % between them and their ideal values is minimized. % Y = INTERP(X,R,L,CUTOFF) allows specification of arguments % L and CUTOFF which otherwise default to 4 and .5 respectively. % 2*L is the number of original sample values used to perform the % interpolation. Ideally L should be less than or equal to 10. % The length of B is 2*L*R+1. The signal is assumed to be band % limited with cutoff frequency 0 < CUTOFF <= 1.0. % [Y,B] = INTERP(X,R,L,CUTOFF) returns the coefficients of the % interpolation filter B. % % % Example: % % Interpolate a signal by a factor of four. % % t = 0:0.001:.029; % Time vector % x = sin(2*pi*30*t) + sin(2*pi*60*t); % Original Signal % y = interp(x,4); % Interpolated Signal % subplot(211); % stem(x); % title('Original Signal'); % subplot(212); % stem(y); % title('Interpolated Signal'); % % See also DECIMATE, RESAMPLE, UPFIRDN. % Author(s): L. Shure, 5-14-87 % L. Shure, 6-1-88, 12-15-88, revised % Copyright 1988-2014 The MathWorks, Inc. % References: % "Programs for Digital Signal Processing", IEEE Press % John Wiley & Sons, 1979, Chap. 8.1. if nargin < 3 l = 4; end if nargin < 4 cutoff = .5; end % Check the input data type. Single precision is not supported. try chkinputdatatype(idata,r,l,cutoff); catch ME throwAsCaller(ME); end if l < 1 || r < 1 || cutoff <= 0 || cutoff > 1 error(message('signal:interp:InvalidRange')) end if abs(r-fix(r)) > eps error(message('signal:interp:MustBeInteger')) end if 2*l+1 > length(idata) s = int2str(2*l+1); error(message('signal:interp:InvalidDimensions', s)); end % ALL occurrences of sin()/() are using the sinc function for the % autocorrelation for the input data. They should all be changed % consistently if they are changed at all. % 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(cutoff*pi*s1)./(cutoff*pi*s1); s2 = sin(cutoff*pi*s2)./(cutoff*pi*s2); ap = s1 + s2; am = s1 - s2; % Compute matrix inverses using Cholesky decomposition for more robustness U = chol(ap); ap = inv(U)*inv(U).'; U = chol(am); am = inv(U)*inv(U).'; % 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:r-1)/r; 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,r) + X*x; y = X + y + eps; h = .5*d'*(sin(pi*cutoff*XX)./(cutoff*pi*XX)); b = zeros(2*l*r+1,1); b(1:l*r) = h'; b(l*r+1) = .5*d(:,l)'*(sin(pi*cutoff*y)./(pi*cutoff*y)); b(l*r+2:2*l*r+1) = b(l*r:-1:1); % use the filter B to perform the interpolation [m,n] = size(idata); nn = max([m n]); if nn == m odata = zeros(r*nn,1); else odata = zeros(1,r*nn); end odata(1:r:nn*r) = idata; % Filter a fabricated section of data first (match initial values and first derivatives by % rotating the first data points by 180 degrees) to get guess of good initial conditions % Filter length is 2*l*r+1 so need that many points; can't duplicate first point or % guarantee a zero slope at beginning of sequence od = zeros(2*l*r,1); od(1:r:(2*l*r)) = 2*idata(1)-idata((2*l+1):-1:2); [od,zi] = filter(b,1,od); %#ok [odata,zf] = filter(b,1,odata,zi); odata(1:(nn-l)*r) = odata(l*r+1:nn*r); % make sure right hand points of data have been correctly interpolated and get rid of % transients by again matching end values and derivatives of the original data if nn == m od = zeros(2*l*r,1); else od = zeros(1,2*l*r); end od(1:r:(2*l)*r) = 2*idata(nn)-(idata((nn-1):-1:(nn-2*l))); od = filter(b,1,od,zf); odata(nn*r-l*r+1:nn*r) = od(1:l*r);