Home > wafo > misc > convlv.m

# convlv

## PURPOSE

Convolves real data set with a response function.

## SYNOPSIS

y = convlv(x,h,dim,flag)

## DESCRIPTION

``` CONVLV Convolves real data set with a response function.

CALL:  y = convlv(x,h,dim,flag);

y,x  = filtered data and raw data, respectively.
h    = vector with filter coefficients.
dim  = dimension to filter.
flag = 'periodic'    : if X is periodic
'nonperiodic' : otherwise (default)

CONVLV convolves X with H, i.e., :

nr
Y(i) = sum Cn(n)*X(i+n)
n=-nl
where
H(1) = Cn(0), H(2)=Cn(-1),...H(Nl+1)=Cn(-Nl), and H(Np) = Cn(1),
H(Np-1) = Cn(2),...H(Np-Nr) = Cn(Nr),

Note: The filter, H, must be stored in wrap-around order, i.e., the
first half of the array H contains the impulse response
function at positive times, while the second half of the array
contains the impulse response function at negative times,
counting down from the highest element.

Example:
dx = 2*pi/100;
x = linspace(0,2*pi-dx,100)';
y = cos(x);
c = savgol(2,2,4,1);          % differentiation filter
yf = convlv(y,c)/dx;          % derivative
yf2 = convlv(y,c,[],'p')/dx;
semilogy(x,abs(sin(x)+yf),x,abs(sin(x)+yf2),'g')

## CROSS-REFERENCE INFORMATION

This function calls:
 conv Convolution and polynomial multiplication. error Display message and abort function. fft Discrete Fourier transform. ifft Inverse discrete Fourier transform. ipermute Inverse permute array dimensions. nan Not-a-Number. nextpow2 Next higher power of 2. permute Permute array dimensions. shiftdim Shift dimensions. strncmpi Compare first N characters of strings ignoring case.
This function is called by:

## SOURCE CODE

```001 function y = convlv(x,h,dim,flag)
002 %CONVLV Convolves real data set with a response function.
003 %
004 %  CALL:  y = convlv(x,h,dim,flag);
005 %
006 %  y,x  = filtered data and raw data, respectively.
007 %  h    = vector with filter coefficients.
008 %  dim  = dimension to filter.
009 %  flag = 'periodic'    : if X is periodic
010 %         'nonperiodic' : otherwise (default)
011 %
012 % CONVLV convolves X with H, i.e., :
013 %
014 %              nr
015 %       Y(i) = sum Cn(n)*X(i+n)
016 %              n=-nl
017 %  where
018 %    H(1) = Cn(0), H(2)=Cn(-1),...H(Nl+1)=Cn(-Nl), and H(Np) = Cn(1),
019 %    H(Np-1) = Cn(2),...H(Np-Nr) = Cn(Nr),
020 %
021 % Note: The filter, H, must be stored in wrap-around order, i.e., the
022 %       first half of the array H contains the impulse response
023 %       function at positive times, while the second half of the array
024 %       contains the impulse response function at negative times,
025 %       counting down from the highest element.
026 %
027 % Example:
028 % dx = 2*pi/100;
029 % x = linspace(0,2*pi-dx,100)';
030 % y = cos(x);
031 % c = savgol(2,2,4,1);          % differentiation filter
032 % yf = convlv(y,c)/dx;          % derivative
033 % yf2 = convlv(y,c,[],'p')/dx;
034 % semilogy(x,abs(sin(x)+yf),x,abs(sin(x)+yf2),'g')
035 %
037
038 %History
039 % Revised pab 2003
040 % - fixed a bug when nr=nl=0 and when nh~=1
041 % By Per A. Brodtkorb 2001
042
043 error(nargchk(2,4,nargin))
044 if nargin<3,dim =[];end
045 if nargin<4|isempty(flag),flag ='nonperiodic';end
046
047 perm = []; nshifts = 0;
048 if ~isempty(dim)
049   perm = [dim:max(ndims(x),dim) 1:dim-1];
050   x = permute(x,perm);
051 else
052   [x,nshifts] = shiftdim(x);
053 end
054 xsiz = size(x);
055
056 if isempty(perm)
057   h    = shiftdim(h);
058   hsiz = size(h);
059 else
060   hsiz = size(h)
061   if prod(hsiz)==prod(xsiz)
062     h = permute(h,perm);
063     hsiz = size(h)
064   end
065 end
066
067 m  = xsiz(1);
068 mh = hsiz(1);
069 nh = prod(hsiz(2:end));
070 p  = ceil(mh/2);
071 nhr = prod(xsiz(2:end)) - nh + 1;
072
073 if (p==mh)
074   y  = x(:,:).*h(ones(m,1),:);
075 else
076   nl = max(find(h(1:p,1)));      % Number of non-zero elements to the left
077   nr = max(find(h(mh:-1:p+1,1)));% and to the right.
078   nz = nr+nl;                    % Total number of non-zero elements
079
080
081   if strncmpi(flag,'p',1)
082     % periodic endconditions
083     nfft = max(m,nz);
084     hw   = fft([h(1:p,:); zeros(nfft-mh,nh); h(p+1:end,:)]);
085     y    = real(ifft(fft(x(:,:),nfft).*repmat(hw,[1,nhr]))); % convolution
086     y    = reshape(y(1:m,:),[ones(1,nshifts),xsiz]);
087   else
088     % nonperiodic end conditions
089     % Extrapolate x linearly outside outside the range to lessen the end
090     % effects
091     y  = linear((1:m)',x(:,:),(-nl+1:m+nr)');
092     %  y  = spline((1:m)',x(:,:),(-nl+1:m+nr)');
093     if 0,
094       %This is wrong
095       y=(conv(y,-h(1:p,:))+flipud(conv(flipud(y),-[0; h(end:-1:p+1,:)])));
096     else
097       nfft = 2^nextpow2(m+nz);
098       hw   = fft([h(1:p,:); zeros(nfft-mh,nh); h(p+1:end,:)]);
099       y  = real(ifft(fft(y(:,:),nfft).*repmat(hw,[1,nhr]))); % convolution
100     end
101     y = reshape(y(nl+1:nl+m,:),[ones(1,nshifts),xsiz]);
102
103   end
104 end
105 if ~isempty(perm), y = ipermute(y,perm); end
106
107
108
109 function yi = linear(x,y,xi);
110 % LINEAR Interpolation and extrapolation
111 %
112 % CALL: yi = linear(x,y,xi);
113 %
114 %  x, xi   = column vectors
115 %  y, yi   = column vectors or matrices
116 %
117
118 siz = size(xi);
119 n = length(x);
120 ni = length(xi);
121 if ni~=1
122     [xxi,k] = sort(xi);
123     [dum,j] = sort([x;xxi]);
124     r(j) = 1:n+ni;
125     r    = r(n+1:end)-(1:ni);
126     r(k) = r;
127     r    = max(1,min(r,n-1));
128     u    = (xi-x(r))./(x(r+1)-x(r));
129     yi   = y(r,:)+(y(r+1,:)-y(r,:)).*u(:,ones(1,size(y,2)));
130 else
131     % Special scalar xi case
132     r = max(find(x <= xi));
133     if isempty(r),  % xi<x(1)
134       r=1;
135     elseif r==n,    % x(n)<=xi
136       r = n-1;
137     end
138     if (r>0) & (r<n),
139         u = (xi-x(r))./(x(r+1)-x(r));
140         yi=y(r,:)+(y(r+1,:)-y(r,:)).*u(:,ones(1,size(y,2)));
141     else
142         yi = NaN;
143     end
144 end
145 if (min(size(yi))==1) & (prod(siz)>1), yi = reshape(yi,siz); end
146
147 return```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003