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') 
  
  See also  savgol

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

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 % 
036 % See also  savgol   
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