Find Optimal Kernel

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

Description

%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

• Required Products: Signal Processing Toolbox
• MATLAB release: MATLAB 7.7 (R2008b)

[Ngawang Jangchub]

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

[Matthew Nelson]

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]

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.

[Ngawang Jangchub]

 

[Ngawang Jangchub]

% 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]

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

[Bart Geurten]

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

[Lakshmi]

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.

[Matthew Nelson]

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.