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
Convolution and polynomial multiplication. | |
Display message and abort function. | |
Discrete Fourier transform. | |
Inverse discrete Fourier transform. | |
Inverse permute array dimensions. | |
Not-a-Number. | |
Next higher power of 2. | |
Permute array dimensions. | |
Shift dimensions. | |
Compare first N characters of strings ignoring case. |
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
Comments or corrections to the WAFO group