Matched Filter

Forum for discussion on different signal processing algorithms
Locked
aloplop
Posts: 41
Joined: 03 Sep 2008, 07:20

Matched Filter

Post by aloplop » 27 Oct 2008, 08:24

Hello,

I wonder if there is any filter for BCI2000 as the one described in the article "A mu-Rhythm Matched Filter for Continuous Control of a Brain-Computer Interface" by Krusienski, Schalk, Mcfarland an Wolpaw. where a matched filter was used for a one-dimensional cursor control task.

Thanks.

mellinger
Posts: 1208
Joined: 12 Feb 2003, 11:06

Post by mellinger » 29 Oct 2008, 07:49

Álvaro,

a matched filter may be implemented using the FIR filter component which is available in the contributions section of BCI2000:
http://www.bci2000.org/wiki/index.php/C ... :FIRFilter

To use it as a replacement for the AR spectral estimator, add FIRFilter.cpp and getfir.cpp to the ARSignalProcessing project in the Borland IDE, and replace the reference to the ARFilter in PipeDefinition.cpp with a reference to the FIRFilter, then recompile and link the module.

To implement a matched filter signal processing, you might want to use spatial filtering in order to extract the source of interest into a single channel, and then use the FIR filter to match it against a mu rhythm pattern. Enter the desired pattern as a single row into the FIRFilterKernal parameter, and set the FIRIntegration parameter to 2 (RMS).

Hope this helps,
Juergen

aloplop
Posts: 41
Joined: 03 Sep 2008, 07:20

Example??

Post by aloplop » 03 Nov 2008, 10:54

Hi,

could you give us an example for that pattern?

In that parameter in the FIRFilter.cpp appears:

Code: Select all

 FIRFilter matrix FIRFilterKernal= 4 4 "
      " 1 0 0 0"
      " 0 1 0 0"
      " 0 0 1 0"
      " 0 0 0 1"
             " 64 -100 100  // Fir Filter Kernal Weights",

So what should we write in the single row for a mu/beta rhythm ??

Thanks.

mellinger
Posts: 1208
Joined: 12 Feb 2003, 11:06

Post by mellinger » 03 Nov 2008, 11:26

Àlvaro,

to extract a mu rhythm wave pattern from an initial recorded session, main analysis steps are:
1) Extract the mu rhythm by an appropriate spatial filter (e.g., Large Laplacian at C3 or C4),
2) Identify the mu rhythm's fundamental frequency,
3) Perform an analysis of relative phase and amplitude between the fundamental frequency and its first two harmonics,
4) Use relative phase and amplitude to construct a wave pattern from sinusoids that correspond to fundamental frequency and harmonics,
5) Enter that wave pattern's sample values into the first row of the FIRFilterKernal parameter.

A quite detailed description of steps 1)-4) is provided in the paper you refer to: Krusienski et al. A mu-Rhythm Matched Filter for Continuous Control of a Brain-Computer Interface. Biomedical Engineering (2007).

Regards,
Juergen
Last edited by mellinger on 04 Nov 2008, 06:05, edited 1 time in total.

aloplop
Posts: 41
Joined: 03 Sep 2008, 07:20

Wave Pattern

Post by aloplop » 03 Nov 2008, 12:04

Hi Juergen,

my colleagues programmed some matlab files to do so for the data in BCI Competition 2003 (data set 3), so I suppose I will have to adapt those files to obtain the wave pattern.

Now my problem is: how to load BCI 2000 data into Matlab and separate between sessions, trials, etc ???

I am using a windows XP 32-bit platform and the Matlab R2007b.

- I don´t know if I can use the mex files --> http://www.bci2000.org/wiki/index.php/U ... _MEX_Files.
I have also downloaded the borland c++ 5.5 compiler (command line version), but typing this doesn´t help:

Code: Select all

make c:\bci2000\tools\mex
- How can I do this if the mex files aren´t adequated ???

Thanks again.

mellinger
Posts: 1208
Joined: 12 Feb 2003, 11:06

Post by mellinger » 04 Nov 2008, 06:04

Hi Álvaro,

in response to your question, I modified the MEX file wiki page, so it should be clearer now.

You don't need to compile MEX files for 32-bit Windows platforms. Just use the binary versions that come with the BCI2000 binary distribution (http://www.bci2000.org/wiki/index.php/BCI2000_Binaries).
Now my problem is: how to load BCI 2000 data into Matlab and separate between sessions, trials, etc ???
A session is a collection of data files that have been recorded in sequence, and are located within a single session directory. Both data file names and directory name reflect the session's number as entered into the SessionNumber parameter. Thus, data is always organized into sessions, and there is no need to separate data into sessions.

A run is the content of a single data file within a session. Thus, data is always separated into runs, and there is no need to separate them.

A trial is a subdivision within a run.

For data recorded using a cursor feedback paradigm, trials are marked by the TargetCode state variable. Each trial corresponds to a period of time where TargetCode is nonzero. Thus, you can identify a start of a new trial by TargetCode changing from zero to nonzero, and the end of a trial by TargetCode changing from nonzero to zero, as in the following Matlab example:

Code: Select all

trial_begin=find((states.TargetCode~=circshift(states.TargetCode,1)).*(states.TargetCode~=0));
trial_end=find((states.TargetCode~=circshift(states.TargetCode,1)).*(states.TargetCode==0));
If you are interested in feedback intervals rather than full trials, simply replace all references to TargetCode with Feedback in the example.

For data recorded using a stimulus presentation paradigm, such as an initial mu rhythm or P300 session, the BCI2000 state variable in question is StimulusCode rather than TargetCode.

Best,
Juergen

aloplop
Posts: 41
Joined: 03 Sep 2008, 07:20

Parameters

Post by aloplop » 11 Nov 2008, 13:49

Hi again,

I have done some analysis of my signals (at the moment I´m not quite sure if they are totally OK). In order to obtain the parameters I do the following:

1º) Load the signals from the Stimulus Presentation. Select the channel where I obtain the best r^2. Do a simple mean between the 5 or 6 runs.

2º) Calculate the Mu rhythm freq.(using the pwelch function). and create a sine at that freq.

Code: Select all

[Pxx fc]=pwelch(means,128,64,256,fs,'onesided');

indice=find((fc>=finfmu)&(fc<fsupmu));

[maxmu indmu]=max(Pxx(indice));       
f_fundmu=fc(min(indice)+indmumu-1);  % Fundamental mu frequency
3º) This sine has 128 samples(= Sampling rate) and I correlate it with portions of the mean signal obtained in (1) with 128 samples. Then I select the portion with most correlation in order to obtain the parameters for the filter.

*** The length of means is a multiple of 128

Code: Select all

for i=1:length(means)/128-1
   aux=means((i-1)*128+1:i*128+1); 
   corrmu=xcorr(aux,mu_sine);
   [m(i) imu(i)]=max(corrmu);
end

In order to do that I normalize this portion and use the fft matlab command with 256 points.

4º) Using the values from the FFT I obtain the amplitudes and phases for the filter.

Code: Select all

TF_norm_means=fft(signal,NFFT)/128;
abs_TF=abs(TF_norm_means);
phase_TF=angle(TF_norm_means);  

% Calculate the parameters of the matched filter
fi1=phase_TF(F==f_fundmu);
fi2=phase_TF(F==2*f_fundmu);
fi3=phase_TF(F==3*f_fundmu);
fi2fi1= fi2-fi1;
fi3fi1= fi3-fi1;

a1=abs_TF(F==f_fundmu);
a2=abs_TF(F==2*f_fundmu);
a3=abs_TF(F==3*f_fundmu);
a2a1=a2/a1;
a3a1=a3/a1;

for n=1:128
 MF(n)=a1*cos(2*pi*(n-1)*1*f_fundmu/fs+fi1)+a2*...
       cos(2*pi*(n-1)*2*f_fundmu/fs+fi2)+a3*cos(2*pi*(n-1)*...
       3*f_fundmu/fs+fi3);
end

MF_norm=MF/max(abs(MF));   % Normalized matched filter.
Is this process done well?? Is there any way to check this?
----------------------------------------------------------------------------------


--> So, if this process were correct what parameters and how should I introduce them in the FIRFIlterkernal parameter to perform the filtering ??


Thanks again.

Alvaro.

mellinger
Posts: 1208
Joined: 12 Feb 2003, 11:06

Post by mellinger » 12 Nov 2008, 06:24

Basically, your method should work. However, there are some disadvantages in the way you implemented it.

You should not average data in step 1), as this will destroy the periodical signal you want to extract. Rather, run the pwelch function over all data, without averaging first. It may also help to constrain this analysis to trials where the mu rhythm is maximally synchronized. As described in the Krusienski paper, averaging should be done after phase alignment (your step 3).

Then, identifying the mu rhythm's fundamental frequency as the spectral maximum within a certain range might lead to false positives, e.g. from noise, muscular artifacts, or a steep slope towards zero. Typically, visual inspection is necessary to determine whether a result makes sense. Fully automatic identification of mu rhythm features is a quite complex problem, and no practically viable method has been established so far.
--> So, if this process were correct what parameters and how should I introduce them in the FIRFIlterkernal parameter to perform the filtering ??
Use the "Normalized matched filter" computed by your last code snippet as a row in the FIRFilter's "FIRFilterKernal" parameter.

HTH,
Juergen

aloplop
Posts: 41
Joined: 03 Sep 2008, 07:20

Post by aloplop » 13 Nov 2008, 13:53

Hello again,

thanks for your recomendations. I have done what you have said and I can see in the plots a signal which actually fits, so I think I am in the right way.

My colleague obtained the best results using a signal segment of 0.25 seconds so in my case I will put in the FIR Windows parameter = 4 (because I use a Sample Block Size = 8 and Sampling rate = 128 Hz).

I would like to ask you some new things:

1º) Applying this filter is like applying a window to the signal so, is the bin width just the same as you told me in http://bci2000.org/phpbb/viewtopic.php?t=418 for the FFTFilter???

I mean, in my case a 32 sample filter corresponds to a 4 Hz width.

2º) What should I do to obtain an adecuate result?? According to the answer of 1º)...

A) Use the 128 sample filter to filter the 32 sample signal.
B) Use the 32 first samples of the filter to compute each 32 block of the signal ??


3º) Is there any manner to pass the coefficients from Matlab to the Matrix in the FIR Filtering tab??

4º) As a small contribution, in the FIRFilter.cpp file the maximum number for integration is 1 but actually should be 3. I had to modify it in order to select the rms integration (value 2).

Thanks.

mellinger
Posts: 1208
Joined: 12 Feb 2003, 11:06

Post by mellinger » 17 Nov 2008, 08:19

ad 1) and 2)
You are right that applying an FIR filter will lead to windowing effects. In your case, the window width needs to be an integer multiple of your signal template, e.g. one to three repetitions of the extracted template.
Also, using just the template as a filter will make the resulting signal jitter because the brain signal will go into phase, out of phase, etc. To avoid this, one would typically use a complex-valued signal template, using exp(iwt) rather than sin(wt) as basis functions. The FIR filter does not allow for complex-valued coefficients, so you might consider building a second template, this time from cos(wt) functions, feed it the same signal as the first one from the spatial filter, and add its output to the first one's in the classifier filter. This way, the BCI control signal should not depend on the arbitrary phase offset between the brain's mu rhythm and your template.

ad 3)
Write the FIRFilterKernal parameter to a file in BCI2000 parameter line format as described in
http://www.bci2000.org/wiki/index.php/T ... eter_Lines
When loading this parameter file fragment from the operator module, only the FIRFilterKernal parameter will be modified.

ad 4)
Thanks for reporting!

aloplop
Posts: 41
Joined: 03 Sep 2008, 07:20

Post by aloplop » 18 Nov 2008, 12:39

Hello,

thanks for your advice...but I am not sure about what I should do.. sorry :(

-> About the filter:

As you have commented, in order to avoid the jitter produced by the phase of the brain signal I should create another "Matched Filter".
In the version I have, cos(wt) functions are used, so I should use the sin() function for this new one and then, normalize it (if I have understood right).

The next I must do is filter the same signal through these two filters.

1º) How can I do this ?? Would be valid if I introduce a second row in the FIRFilteKernal parameter with this new matched filter??

Let´s say, for example, I have 64 samples of the brain signal and, although both matched filters have 128 samples (as I use a sampling freq. of 128 Hz), I use the first 64 samples of both filters in the FIR Filter:

2º) How can I know what samples should I introduce in the Classifier???
This aspect depends on how the previous computation is performed (i.e. how the signal is filtered by the FIR).


--> About the parameter line:

I have created a .prm file and it works properly just by loading the file created by you.

Thanks again.

mellinger
Posts: 1208
Joined: 12 Feb 2003, 11:06

Post by mellinger » 19 Nov 2008, 07:24

ad 0)
exp(iwt) = cos(wt) + i*sin(wt)
Thus, to get real and imaginary part as two different filters, use cosines for one, and sines for the other. If you would like to normalize them (I don't see why you would bother), both should be multiplied using the same normalization factor.

ad 1)
Yes, introduce a second row for the second filter, and make sure there exists an input channel for this row, by duplicating the first row of your spatial filter matrix into the second one. Again, the length of the FIR filter should match the length of your template, or an integer multiple thereof, not the length of a sampling block. In case the FIRFilter doesn't allow for this, simply pad with zeros at the end.

ad 2)
Try the first sample in the block.

HTH,
Juergen

aloplop
Posts: 41
Joined: 03 Sep 2008, 07:20

Post by aloplop » 10 Dec 2008, 06:32

Hello,

Finally, I decided to use a Matched FIR Filter of 32 samples.

However, when I have to correlate a piece of the signal with a sine at the mu rhythm fundamental frequency I don´t know what would be the right length (in samples) of this sine in order to obtain the pattern for the matched filters.

I don´t know whether to take 32 (length for the filter) or 128 samples (sampling frequency).

I have performed a test where it is better to take 128 samples for the sine. Maybe it is because " the mu rhythm tends to appear in bursts lasting anywhere from less than 1s to several seconds", as described in the article of Krusienski et. al.

Thanks again.

mellinger
Posts: 1208
Joined: 12 Feb 2003, 11:06

Post by mellinger » 10 Dec 2008, 08:27

Álvaro,

as an example, consider a sampling frequency of 128 Hz, and a fundamental mu rhythm frequency of 10 Hz. Then, a full period of the sine will correspond to 12.8 samples, and your pattern should have a length of 13, 26, 38, ... samples.

Using a number of repeated waveforms is equivalent to running a boxcar-like lowpass filter over the FIRFilter's output.
You might consider using a single waveform only, and do subsequent low-pass filtering by means of the LPFilter, or alternatively a number of waveforms, and not use the LPFilter. When using a number of waveform repetitions, it might also make sense to weight FIR coefficients with a window (Hann, Hanning, Blackman ...) as in case of the FFT, to reduce windowing effects (but I doubt that it will make much of a difference).

Without knowing your goodness criterion, I cannot comment on the result that "128 samples is better than 32". At a sampling rate of 128 Hz, 32 samples correspond to a frequency of 4 Hz, so I don't think this actually matches the mu rhythm's fundamental frequency, and consequently I wouldn't expect it to work well.

HTH,
Juergen

Locked

Who is online

Users browsing this forum: No registered users and 4 guests