nt_sparse_filter

y=nt_sparse_filter(x,T,A) - convolve multichannel data with sparse impulse response

  y: result

  x: data to convolve (columnwise)
  T: times of non-zero samples of IR (can be fractionnary)
  A: amplitude of non-zero samples of IR (default: all 1)

 T and A together describe the impulse response.  A must have same size
 as T.

 If T is a column vector it is applied to all columns of x.

 If T is a 2D matrix or 1D cell array each column or cell is applied to a
 column of x (dimensions must fit).

 If T is a 2D cell array it is interpreted as defining a multichannel 
 impulse response. The cell on row I and column J is applied to column I
 of x and the result is added to column J of y.

 Fractionary lags are implemented by linear interpolation.

0001 function y=nt_sparse_filter(x,T,A)
0002 %y=nt_sparse_filter(x,T,A) - convolve multichannel data with sparse impulse response
0003 %
0004 %  y: result
0005 %
0006 %  x: data to convolve (columnwise)
0007 %  T: times of non-zero samples of IR (can be fractionnary)
0008 %  A: amplitude of non-zero samples of IR (default: all 1)
0009 %
0010 % T and A together describe the impulse response.  A must have same size
0011 % as T.
0012 %
0013 % If T is a column vector it is applied to all columns of x.
0014 %
0015 % If T is a 2D matrix or 1D cell array each column or cell is applied to a
0016 % column of x (dimensions must fit).
0017 %
0018 % If T is a 2D cell array it is interpreted as defining a multichannel
0019 % impulse response. The cell on row I and column J is applied to column I
0020 % of x and the result is added to column J of y.
0021 %
0022 % Fractionary lags are implemented by linear interpolation.
0023 
0024 if nargin<3; A=[]; end
0025 if nargin<2; error('!'); end
0026 
0027 % check A or create
0028 %if isvector(T); T=T(:); A=A(:); end
0029 if isnumeric(T)
0030     if isempty(A); A=ones(size(T)); end
0031     if ~isnumeric(A); error('!'); end
0032     if size(A) ~= size(T); error('!'); end
0033     if size(A,2) ~= 1 || (size(A,2) ~= size(x,2) && size(A,2)~=1); error('!'); end
0034 elseif iscell(T)
0035     if isempty(A); 
0036         A=cell(size(T));
0037         for iCell=1:numel(T); A{iCell}=ones(size(T{iCell})); end
0038     end
0039     if ~iscell(A); error('!'); end
0040     if size(A) ~= size(T); error('!'); end
0041     if size(A,1)~=size(x,2); error('number of rows of cell array of IRs should equal number of columns of x'); end
0042 end
0043 
0044 
0045 % handle negative lags
0046 if isnumeric(T)
0047     a=min(T(:));
0048     if a<0
0049         x=[zeros(ceil(-a)),size(x,2)];
0050         T=T+ceil(-a);
0051     end
0052 elseif iscell(T)
0053     a=min(T{1,1}); 
0054     for iCell=1:numel(T)
0055         a=min(a,min(T{iCell}));
0056     end
0057     if a<0
0058         x=[zeros(ceil(-a)),size(x,2)];
0059         for iCell=1:numel(T)
0060             T{iCell}=T{iCell}+ceil(-a);
0061         end
0062     end
0063 else; error('!'); end % neither matrix nor cell
0064 
0065 % estimate nrows of result
0066 if isnumeric(T)
0067     b=max(T(:));
0068 else
0069     b=max(T{1,1});
0070     for iCell=1:numel(T)
0071         b=max(b,max(T{iCell}));
0072     end
0073 end
0074 b=ceil(b);
0075 
0076 [nsamples,ncolsx]=size(x);
0077 
0078 if isnumeric(T) && size(T,2)==1 
0079     % apply same IR to all columns of x
0080     y=zeros(nsamples+b, ncolsx);
0081     for iPulse=1:numel(T);
0082         t=T(iPulse);
0083         integT=floor(t); 
0084         fracT=t-integT;
0085         y(integT+(1:nsamples),:) = y(integT+(1:nsamples),:) + x*(1-fracT)*A(iPulse);
0086         if fracT; y(1+integT+(1:nsamples),:) = y(1+integT+(1:nsamples),:) + x*fracT*A(iPulse); end
0087     end
0088 end
0089 
0090 if isnumeric(T) && size(T,2)>1
0091     % apply different IR to each column of x
0092     y=zeros(nsamples+b, ncolsx);
0093     for iCol=1:ncolsx
0094         for iPulse=1:numel(T(:,iCol));
0095             t=T(iPulse,iCol);
0096             integT=floor(t); 
0097             fracT=t-integT;
0098             1-fracT
0099             y(integT+(1:nsamples),iCol) = y(integT+(1:nsamples),iCol) + x(:,iCol)*(1-fracT)*A(iPulse,iCol);
0100             if fracT; y(1+integT+(1:nsamples),iCol) = y(1+integT+(1:nsamples),iCol) + x(:,iCol)*fracT*A(iPulse,iCol); end
0101         end
0102     end
0103 end
0104 
0105 if iscell(T) && numel(T)==ncolsx
0106     % apply one IR to each column, no cross-terms
0107     y=zeros(nsamples+b, ncolsx);
0108     for iCol=1:ncolsx
0109         for iPulse=1:numel(T{iCol});
0110             t=T{iCol}(iPulse);
0111             integT=floor(t); 
0112             fracT=t-integT;
0113             y(integT+(1:nsamples),iCol) = y(integT+(1:nsamples),iCol) + x(:,iCol)*(1-fracT)*A{iCol}(iPulse);
0114             if fracT; y(1+integT+(1:nsamples),iCol) = y(1+integT+(1:nsamples),iCol) + x(:,iCol)*fracT*A{iCol}(iPulse); end
0115         end
0116     end
0117 end
0118 
0119 if iscell(T) && size(T,1)==ncolsx 
0120     ncolsy=size(T,2);
0121     % full multichannel IR with cross-terms
0122     y=zeros(nsamples+b, ncolsy);
0123     for iRow=1:ncolsx
0124         for iCol=1:ncolsy
0125             for iPulse=1:numel(T{iRow,iCol})
0126                 t=T{iRow,iCol}(iPulse);
0127                 integT=floor(t); 
0128                 fracT=t-integT;
0129                 y(integT+(1:nsamples),iCol) = y(integT+(1:nsamples),iCol) + x(:,iRow)*(1-fracT)*A{iRow,iCol}(iPulse);
0130                 if fracT; y(1+integT+(1:nsamples),iCol) = y(1+integT+(1:nsamples),iCol) + x(iRow)*fracT*A{iRow,iCol}(iPulse); end
0131             end
0132         end
0133     end
0134 end
0135 
0136 
0137 return
0138 
0139 % test code
0140 x=randn(1000,1);
0141 y=nt_sparse_filter(x,[0]);
0142 figure(1); clf; subplot 211; plot(x); subplot 212; plot(y);
0143 
0144 x=randn(1000,1);
0145 y=nt_sparse_filter(x,[0.5]);
0146 figure(1); clf; subplot 211; plot(x); subplot 212; plot(y);
0147 
0148 x=randn(1000,1);
0149 y=nt_sparse_filter(x,[10]);
0150 figure(1); clf; subplot 211; plot(x); subplot 212; plot(y);
0151 
0152 x=randn(1000,1);
0153 y=nt_sparse_filter(x,[1:10]);
0154 figure(1); clf; subplot 211; plot(x); subplot 212; plot(y);
0155 
0156 x=randn(1000,1);
0157 y=nt_sparse_filter(x,[1:10], ones(1,10));
0158 figure(1); clf; subplot 211; plot(x); subplot 212; plot(y);
0159 
0160 x=randn(1000,2);
0161 y=nt_sparse_filter(x,[1:10], ones(1,10));
0162 figure(1); clf; subplot 211; plot(x); subplot 212; plot(y);
0163 
0164 x=randn(1000,2);
0165 y=nt_sparse_filter(x,[1:10], ones(2,10));
0166 figure(1); clf; subplot 211; plot(x); subplot 212; plot(y);
0167 
0168 x=randn(1000,2);
0169 y=nt_sparse_filter(x,{ones(1,10),ones(1,10)});
0170 figure(1); clf; subplot 211; plot(x); subplot 212; plot(y);
0171 
0172 x=randn(1000,1);
0173 y=nt_sparse_filter(x,{ones(1,10),ones(1,10)});
0174 figure(1); clf; subplot 211; plot(x); subplot 212; plot(y);
0175 
0176             
0177