Home > NoiseTools > nt_sparse_filter.m

# nt_sparse_filter

## PURPOSE

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

## SYNOPSIS

function y=nt_sparse_filter(x,T,A)

## DESCRIPTION

    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.

## CROSS-REFERENCE INFORMATION

This function calls:

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

This function is called by:

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

## SOURCE CODE

    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