File Exchange

## Find Optimal Kernel

version 1.2 (1.76 KB) by

Find optimal causal filter kernel that convolves with one time series to approximate another

Updated

%function k=FindOptKern(InSig,OutSig,MaxnLags)
%
% Given InSig, OutSig, and the MaxnLags, this will find the optimal causal
% kernal of length MaxnLags that when convolved with InSig gives the
% optimal approximation to OutSig in a least squares sense.
%
% Note that this code produces the optimal causal kernel. the basic
% technique does not have to produce a causal kernel, but this code does.
%
% Input
% InSig: An N x 1 vector of the input signal (i.e. Trial type in the
% trial history example)
% OutSig: An N x 1 vector of the output signal (i.e. RT in the trial
% history example)
% MaxnLags: A scalar for the maximum number of lags to calculate up to.
% Must be less than the length of the signal.
%
% k: A MaxnLags x 1 vector of the optimal kernel.
%
% For an explanation of this, See Corrado, Sugrue, Seung and Newsome (2004)
% pp. 591-592
% Also see: http://hebb.mit.edu/courses/9.29/2004/lectures/lecture03.pdf
%
% by MJN on 090318

Luke Gahan

Luke Gahan

### Luke Gahan (view profile)

Hi,

Just one small question. For xcov(x,y) you retain Cxy(MaxnLags+2:end). This means that you are discarding the covariance at zero lag and retaining 1:MaxnLags.

While for xcov(x,x) you retain 0:MaxnLags-1.

My understanding of the procedure (which could well be incorrect) is that if the user specifies MaxnLags, then maybe 0:MaxnLags should be retained, meaning MaxnLags+1 filter coefficients will be calculated (zero time plus lags)?

Thanks,

Luke

Matthew Nelson

### Matthew Nelson (view profile)

These are just parameters in creating a sample kernel for testing the performance of the function. You can adjust them however you please to change the test kernel. The purpose of including TestFindOptKern.m was just to demonstrate the use of the function, and is not by itself an important script.

Lakshmi

### Lakshmi (view profile)

Tried to run the test script that you have provided. But there are some undefined variables like BW, Growth, Decay, etc. Can you please update the file or atleast provide instructions for initializing the undefined variables in the test script.

Bart Geurten

### Bart Geurten (view profile)

I just used your code and it is good commented. Does what I wanted it to do ;)

Ngawang Jangchub

### Ngawang Jangchub (view profile)

% Sorry for the typos above...

newInSig=[InSig; InSig];
newOutSig=[OutSig; OutSig];

for i = 1:r
Cxy(i,1)=InSig' * newOutSig(i+r-1,:);
Cxx(i,1)=InSig' * newInSig(i+r-1,:);
end
Cxy=Cxy(2:MaxnLags+1, :);
Cxx=Cxx(1:MaxnLags, :);
Cxxmat=toeplitz(Cxx);

k= (Cxxmat^-1)*Cxy;

Ngawang Jangchub

### Ngawang Jangchub (view profile)

% Without access to the Signal Processing Toolbox, it's difficult to % say for sure, but I suspect the following code would allow you to % run this routine without calling the xcov() function.

newInSig=[InSig; InSig];
newOutSig=[OutSig;outSig];

for i = 1:r
Cxy(i,1)=InSig' * newOutSig(i+r-1,:);
Cxx(i,1)=InSig' * newOutSig(i+r-1,:);
end
Cxy=Cxy(2:MaxnLags+1, :);
Cxx=Cxx(1:MaxnLags, :);
Cxxmat=toeplitz(Cxx);

k= (Cxxmat^-1)*Cxy;

Ngawang Jangchub

Ngawang Jangchub

### Ngawang Jangchub (view profile)

Maybe you could provide a simple sample input & output so users can know they have the code working properly -- especially those of us who don't have the signal processing toolbox. Thx, N.

Matthew Nelson

### Matthew Nelson (view profile)

xcov is available in matlab's signal processing toolbox. I'm sorry if that wasn't clear before. I've adjusted the online page to reflect that.

Ngawang Jangchub

### Ngawang Jangchub (view profile)

Unless I missed something, you didn't provide the code for the xcov() function. You call this function in FindOptKern()...
Cxy=xcov(OutSig,InSig,MaxnLags);