Home > NoiseTools > nt_mchanfilt.m

# nt_mchanfilt

## PURPOSE

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

## SYNOPSIS

function y=nt_mfilt(x,M,B,A,expand)

## DESCRIPTION

    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);

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    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