y=nt_mfilt(x,M,B,A,expand) - multichannel filter y: filtered data x: data to filter (time X channel) M: multichannel impulse response 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 band, the spatial filter is defined by one page of M, and the temporal filter by one column of B (and A if IIR). Default filters are 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 ) Examples: Filter with multichannel FIR filter M: y=nt_mfilt(x,M) Same, but lags replaced by FIR filters: y=nt_mfilt(x,M,B) Same, but filters are IIR: y=nt_mfilt(x,M,B,A); Examples of filter bases: Basis of lags (default): B=eye(nbands); Basis of nbands cosines of duration 100 samples: B=cos(2*pi*(1:100)'*(1:nbands)/100) Basis of 6 dyadic filters: b=zeros(32,1); B=nt_multismooth(b,[1 2 4 8 16 32],[],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 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 band, the spatial filter is defined by one page of M, and the 0013 % temporal filter by one column of B (and A if IIR). 0014 % 0015 % Default filters are B=eye(nbands) and A=ones(1,nbands); 0016 % 0017 % Data can be 2D matrix or cell array of 2D matrices. 0018 % 0019 % M is 3D matrix ( inchannels X outchannels X bands ) 0020 % 0021 % Examples: 0022 % Filter with multichannel FIR filter M: 0023 % y=nt_mfilt(x,M) 0024 % 0025 % Same, but lags replaced by FIR filters: 0026 % y=nt_mfilt(x,M,B) 0027 % 0028 % Same, but filters are IIR: 0029 % y=nt_mfilt(x,M,B,A); 0030 % 0031 % Examples of filter bases: 0032 % Basis of lags (default): 0033 % B=eye(nbands); 0034 % Basis of nbands cosines of duration 100 samples: 0035 % B=cos(2*pi*(1:100)'*(1:nbands)/100) 0036 % Basis of 6 dyadic filters: 0037 % b=zeros(32,1); B=nt_multismooth(b,[1 2 4 8 16 32],[],1); 0038 % 0039 0040 0041 if nargin<2; error('!'); end 0042 if nargin<3; B=[]; end 0043 if nargin<4; A=[]; end 0044 if nargin<5; expand=0; end 0045 0046 % handle cell array data 0047 if iscell(x) 0048 y={}; 0049 for iCell=1:numel(x) 0050 y{iCell}=nt_mfilt(x{iCell},M,B,A,expand); 0051 end 0052 return; 0053 end 0054 if ndims(x)>2; error('2D only - pass higher dim data as cell array'); end 0055 0056 % sizes consistent? 0057 [nchans_i,nchans_o,nbands]=size(M); 0058 if size(x,2) ~= nchans_i ; error('!'); end 0059 0060 % default filters 0061 if isempty(B); B=eye(nbands); end 0062 if isempty(A); A=ones(1,nbands); end 0063 % check sizes 0064 if size(B,2) ~= nbands; error('!'); end 0065 if size(A,2) ~= nbands; error('!'); end 0066 0067 % filter 0068 y=zeros(size(x,1),nchans_o,nbands); 0069 for iBand=1:nbands 0070 xx=filter(B(:,iBand),A(:,iBand),x); 0071 y(:,:,iBand)=xx*M(:,:,iBand); 0072 end 0073 0074 if ~expand 0075 y=sum(y,3); 0076 end 0077 0078 % tests/examples 0079 if 0 0080 % basic tests 0081 x=rand(100,1); % single channel data 0082 M=ones(1,1,1); 0083 y=nt_mfilt(x,M); 0084 0085 B=1; 0086 y=nt_mfilt(x,M,B); 0087 0088 A=1; 0089 y=nt_mfilt(x,M,B,A); 0090 0091 M=ones(1,1,10); % 10-tap FIR 0092 y=nt_mfilt(x,M); 0093 0094 M=ones(1,5,1); % fanout to 5 channels 0095 y=nt_mfilt(x,M); 0096 0097 M=ones(1,5,10); % fanout to 5, 10-tap FIR 0098 y=nt_mfilt(x,M); 0099 0100 x=randn(100,15); % 15-channel data 0101 M=ones(15,5,1); % fanin to 5 0102 y=nt_mfilt(x,M); 0103 0104 M=ones(15,5,10); % fanin to 5, 10-tap FIR 0105 y=nt_mfilt(x,M); 0106 0107 B=eye(10); % basis is lags (default) 0108 y=nt_mfilt(x,M,B); 0109 0110 B=ones(11,10); % basis is 10-channel filterbank made of FIRs of order 11 0111 y=nt_mfilt(x,M,B); 0112 0113 B=ones(3,10); A=ones(2,10); % basis is 10-channel filterbank made of IIRs of order 3 0114 y=nt_mfilt(x,M,B,A); 0115 end 0116 0117 if 0 0118 x=zeros(100,1); x(1)=1.1; % 0119 M=zeros(1,1,6); M(1,1,6)=1; % delay by 5 samples 0120 figure(1); clf; plot(nt_mfilt(x,M)); 0121 0122 M=zeros(1,6,6); 0123 for k=1:6;M(1,k,k)=1; end; % delay by 0:5 samples 0124 figure(1); clf; plot(nt_mfilt(x,M)); 0125 0126 B=zeros(61,6); 0127 for k=1:6; B((k-1)*10+1,k)=1; end; % basis consists of set of larger delays 0128 figure(1); clf; plot(nt_mfilt(x,M,B)); 0129 0130 B=[]; A=[]; 0131 for k=1:6; 0132 [B(k,:),A(k,:)]=butter(2,[k,k+1]/(2*10),'bandpass'); % basis consists of bandpass filters 0133 end 0134 B=B'; A=A'; 0135 figure(1); clf; plot(nt_mfilt(x,M,B,A)); 0136 0137 x=randn(100,3); % 3-channel 0138 M=randn(3,4,6); % fanout to 4, order-6 'FIR' 0139 y=nt_mfilt(x,M,B,A); % apply using bandpass basis 0140 figure(1); clf; plot(y); 0141 % The output here is the sum of 6 4-channel signals, each produced by 0142 % applying a 3X4 transform matrix to input signals filtered by the 0143 % corresponding basis. 0144 0145 end 0146 0147 0148 if 0 0149 % equivalent of nt_multishift 0150 x=zeros(100,1); 0151 x(1,:)=1; 0152 expand=true; 0153 y=nt_mfilt(x,ones(1,1,10),eye(10),[],expand); 0154 disp(size(y)) 0155 figure(1); clf; plot(squeeze(y*1.1)); 0156 end 0157 0158