nt_mfilt

y=nt_mfilt(x,M,B,A,expand) - multichannel filter

  y: filtered data 

  x: data to filter (time X channel) 
  M: multichannel impulse response (ichannel X ochannel X lag)
  B, A: bank of filters 
  expand: if true output 3D matrix, one band per page (else add)

 Output is sum of spatially and temporally filtered inputs, one for each band. 
 For each lag, the spatial filter is defined by one page of M.

 If B is provided, lags are replaced by FIR filter outputs (each
 column of A defines a FIR filter).  If A is provided, the
 filters are IIR.
 
 Default filters are lags: B=eye(nbands) and A=ones(1,nbands);

  Data can be 2D matrix or cell array of 2D matrices.

  M is 3D matrix ( inchannels X outchannels X bands ). If M is empty, it
  is replaced as:      
       M=zeros(size(x,2),size(x,2),size(B,2));
       for k=1:size(x,2); M(k,k,:)=1;

0001 function y=nt_mfilt(x,M,B,A,expand)
0002 %y=nt_mfilt(x,M,B,A,expand) - multichannel filter
0003 %
0004 %  y: filtered data
0005 %
0006 %  x: data to filter (time X channel)
0007 %  M: multichannel impulse response (ichannel X ochannel X lag)
0008 %  B, A: bank of filters
0009 %  expand: if true output 3D matrix, one band per page (else add)
0010 %
0011 % Output is sum of spatially and temporally filtered inputs, one for each band.
0012 % For each lag, the spatial filter is defined by one page of M.
0013 %
0014 % If B is provided, lags are replaced by FIR filter outputs (each
0015 % column of A defines a FIR filter).  If A is provided, the
0016 % filters are IIR.
0017 %
0018 % Default filters are lags: B=eye(nbands) and A=ones(1,nbands);
0019 %
0020 %  Data can be 2D matrix or cell array of 2D matrices.
0021 %
0022 %  M is 3D matrix ( inchannels X outchannels X bands ). If M is empty, it
0023 %  is replaced as:
0024 %       M=zeros(size(x,2),size(x,2),size(B,2));
0025 %       for k=1:size(x,2); M(k,k,:)=1;
0026 
0027 %
0028 %  Examples:
0029 %  Filter with multichannel FIR filter M:
0030 %       y=nt_mfilt(x,M)
0031 %
0032 %  Same, but lags replaced by FIR filters:
0033 %       y=nt_mfilt(x,M,B)
0034 %
0035 %  Same, but filters are IIR:
0036 %       y=nt_mfilt(x,M,B,A);
0037 %
0038 %  Examples of filter bases:
0039 %  Basis of lags (default):
0040 %       B=eye(nbands);
0041 %  Basis of nbands cosines of duration 100 samples:
0042 %       B=cos(2*pi*(1:100)'*(1:nbands)/100)
0043 %  Basis of 6 dyadic filters:
0044 %       b=zeros(32,1); B=nt_multismooth(b,[1 2 4 8 16 32],[],1);
0045 %
0046 % Other simple examples:
0047 %
0048 %  Simple matrix multiplication:
0049 %       M=zeros(size(M0,1),size(M0,2),1); % M0 is the matrix
0050 %       M(:,:,1)=M0;
0051 %       y=nt_mfilt(x,M);
0052 %
0053 %  Apply a different FIR to each data channel:
0054 %       M=zeros(size(x,2),size(x,2),order);
0055 %       for k=1:size(x,2); M(k,k,:)=B(:,k); % FIRs are columns of B
0056 %       y=nt_mfilt(x,M,B);
0057 %
0058 %  Same, different approach:
0059 %       M=zeros(size(x,2),size(x,2),size(x,2));
0060 %       for k=1:size(x,2); M(k,k,k)=1; end
0061 %       y=nt_mfilt(x,M,B);
0062 %
0063 %  Apply the same set of FIR filters B to each channel
0064 %       M=[];
0065 %       expandflag=1; % output each filter on different page
0066 %       y=nt_mfilt(x,M,B,[],expandflag); % y is 3D
0067 
0068 if nargin<2; error('!'); end
0069 if nargin<3; B=[]; end
0070 if nargin<4; A=[]; end
0071 if nargin<5; expand=0; end
0072 
0073 % handle cell array data
0074 if iscell(x)
0075     y={};
0076     for iCell=1:numel(x)
0077         y{iCell}=nt_mfilt(x{iCell},M,B,A,expand);
0078     end
0079     return;
0080 end
0081 if ndims(x)>2; error('2D only - pass higher dim data as cell array'); end
0082 
0083 if isempty(M)
0084     M=zeros(size(x,2),size(x,2),size(B,2));
0085     for k=1:size(x,2); M(k,k,:)=1; end
0086 end
0087 
0088 % sizes consistent?
0089 [nchans_i,nchans_o,nbands]=size(M);
0090 if size(x,2) ~= nchans_i ; error('!'); end
0091 
0092 % default filters
0093 if isempty(B); B=eye(nbands); end
0094 if isempty(A); A=ones(1,nbands); end
0095 % check sizes
0096 if size(B,2) ~= nbands;  error('!'); end
0097 if size(A,2) ~= nbands;  error('!'); end
0098     
0099 % filter
0100 y=zeros(size(x,1),nchans_o,nbands);
0101 for iBand=1:nbands
0102     xx=filter(B(:,iBand),A(:,iBand),x);
0103     y(:,:,iBand)=xx*M(:,:,iBand);
0104 end
0105 
0106 if ~expand
0107     y=sum(y,3);
0108 end
0109 
0110 % tests/examples
0111 if 0
0112     % basic tests
0113     x=rand(100,1); % single channel data
0114     M=ones(1,1,1);
0115     y=nt_mfilt(x,M);
0116     disp(size(y));
0117     
0118     B=1;
0119     y=nt_mfilt(x,M,B);
0120     disp(size(y));
0121     
0122     A=1;
0123     y=nt_mfilt(x,M,B,A);
0124     disp(size(y));
0125     
0126     M=ones(1,1,10); % 10-tap FIR
0127     y=nt_mfilt(x,M); 
0128     disp(size(y));
0129     
0130     M=ones(1,5,1); % fanout to 5 channels
0131     y=nt_mfilt(x,M); 
0132     disp(size(y));
0133     
0134     M=ones(1,5,10); % fanout to 5, 10-tap FIR
0135     y=nt_mfilt(x,M); 
0136     disp(size(y));
0137     
0138     x=randn(100,15); % 15-channel data
0139     M=ones(15,5,1); % fanin to 5
0140     y=nt_mfilt(x,M); 
0141     disp(size(y));
0142     
0143     M=ones(15,5,10); % fanin to 5, 10-tap FIR
0144     y=nt_mfilt(x,M); 
0145     disp(size(y));
0146     
0147     B=eye(10); % basis is lags (default)
0148     y=nt_mfilt(x,M,B); 
0149     disp(size(y));
0150 
0151     B=ones(11,10); % basis is 10-channel filterbank made of FIRs of order 11
0152     y=nt_mfilt(x,M,B); 
0153     disp(size(y));
0154 
0155     B=ones(3,10); A=ones(2,10); % basis is 10-channel filterbank made of IIRs of order 3
0156     y=nt_mfilt(x,M,B,A);    
0157     disp(size(y));
0158 end
0159     
0160 if 0
0161     x=zeros(100,1); x(1)=1.1; %
0162     M=zeros(1,1,6); M(1,1,6)=1; % delay by 5 samples
0163     figure(1); clf; plot(nt_mfilt(x,M));
0164     
0165     M=zeros(1,6,6); 
0166     for k=1:6;M(1,k,k)=1; end; % delay by 0:5 samples
0167     figure(1); clf; plot(nt_mfilt(x,M));
0168 
0169     B=zeros(61,6); 
0170     for k=1:6; B((k-1)*10+1,k)=1; end;  % basis consists of set of larger delays
0171     figure(1); clf; plot(nt_mfilt(x,M,B));    
0172     
0173     B=[]; A=[];
0174     for k=1:6;
0175         [B(k,:),A(k,:)]=butter(2,[k,k+1]/(2*10),'bandpass'); % basis consists of bandpass filters
0176     end
0177     B=B'; A=A';
0178     figure(1); clf; plot(nt_mfilt(x,M,B,A));   
0179     
0180     x=randn(100,3); % 3-channel
0181     M=randn(3,4,6); % fanout to 4, order-6 'FIR'
0182     y=nt_mfilt(x,M,B,A);  % apply using bandpass basis
0183     figure(1); clf; plot(y);
0184     % The output here is the sum of 6 4-channel signals, each produced by
0185     % applying a 3X4 transform matrix to input signals filtered by the
0186     % corresponding basis.
0187     
0188 end
0189 
0190 
0191 if 0
0192     % equivalent of nt_multishift
0193     x=zeros(100,1);
0194     x(1,:)=1;
0195     expand=true;
0196     y=nt_mfilt(x,ones(2,1,10),eye(10),[],expand);
0197     disp(size(y))
0198     figure(1); clf; plot(squeeze(y*1.1));
0199 end
0200 
0201