Code covered by the BSD License

### Highlights from Find Optimal Kernel

4.0
4.0 | 5 ratings Rate this file 8 Downloads (last 30 days) File Size: 2.63 KB File ID: #23363 Version: 1.2

# Find Optimal Kernel

### Matthew Nelson (view profile)

19 Mar 2009 (Updated )

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

File Information
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)
16 May 2013 Luke Gahan

### Luke Gahan (view profile)

30 Mar 2013 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

Comment only
20 Jan 2012 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.

Comment only
19 Jan 2012 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.

03 Jun 2009 Bart Geurten

### Bart Geurten (view profile)

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

26 Mar 2009 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;

Comment only
25 Mar 2009 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;

Comment only
25 Mar 2009 Ngawang Jangchub

### Ngawang Jangchub (view profile)

25 Mar 2009 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.

Comment only
25 Mar 2009 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.

Comment only
25 Mar 2009 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);