Skip to content

Commit e287299

Browse files
chkothechkothe
committed
Merge pull request #16 from nedaeh/EmotionRecognition2
Modified version of FBCSP using welch method, this replicates the cur…
2 parents 5b7db7f + 27cd14d commit e287299

File tree

1 file changed

+243
-0
lines changed

1 file changed

+243
-0
lines changed

‎code/paradigms/ParadigmFBCSPMod.m

Lines changed: 243 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,243 @@
1+
classdef ParadigmFBCSPMod < ParadigmDataflowSimplified
2+
% Paradigm for complex oscillatory processes using the filter-bank CSP algorithm.
3+
% Result = para_multiband_csp(Input-Data, Operation-Mode, Options...)
4+
%
5+
% Filter-bank CSP [1,2] is a simple extension of the basic CSP method (see ParadigmCSP), in which for
6+
% each of several time/frequency filters a set of CSP filters is learned, followed by log-variance
7+
% feature extraction, concatenation of all features (over all chosen spectral filters) and
8+
% subsequent machine learning. It is not a general replacement for CSP due to the problem of
9+
% overfitting, but is very useful whenever oscillatory processes in different frequency bands (and
10+
% with different spatial topographies) are jointly active, and their concerted behavior must be
11+
% taken into account for a given prediction task. Filter-bank CSP can also be used to capture
12+
% oscillations in multiple time windows, instead of frequency windows (for example for the detection
13+
% of complex event-related dynamics).
14+
%
15+
% Since the dimensionality of the feature space is larger than in CSP, and since complex
16+
% interactions may be present, a more complex classifier than the default LDA may be necessary to
17+
% learn an appropriate model. On the other hand, more flexibility amplifies the risk of overfitting
18+
% (especially with only little calibration data), so that the performance should always be compared
19+
% to standard CSP (and Spec-CSP). Another reason is that complex (relevant) interactions between
20+
% different frequency bands are seemingly rarely observed in practice. The most important
21+
% user-configurable parameters are the selection regions in time and frequency and the learner
22+
% component.
23+
%
24+
% Typical applications would be those in which either complex event-related oscillatory dynamics
25+
% happen (for example when reacting to a particular stimulus) and/or where non-trivial interactions
26+
% between frequency bands (e.g. alpha/theta) are relevant, such as, for example, in workload
27+
% measurements.
28+
%
29+
% Example: Consider a calibration data set in which a subject is maintaining and updating
30+
% different number of items in his/her working memory at different times, e.g. while performing
31+
% the n-Back task [2]. Events with types 'n1','n2','n3' indicate challenge stimuli in which the
32+
% respective number of items is being processed by the person. The goal is to be able to predict
33+
% the working-memory load of the person following the presentation of such a memory-related
34+
% challenge. An epoch of 3 seconds relative to each challenge is selected, and three different
35+
% regions are chosen, two of them over the entire interval, covering the theta and alpha ryhthm,
36+
% respectively, and one region that is restricted to a window around the time of heaviest
37+
% cognitive processing. The three regions are specified as a cell array of flt_select
38+
% parameters.
39+
%
40+
% data = io_loadset('data sets/mary/nback.eeg')
41+
% myapproach = {'FBCSP' 'SignalProcessing',{'EpochExtraction',[-0.5 2.5]}, ...
42+
% 'Prediction', {'FeatureExtraction',{'FreqWindows',[4 6; 7 15; 7 15],'TimeWindows',[-0.5 2.5; -0.5 2.5; 0.25 1.25]}, ...
43+
% 'MachineLearning',{'Learner','logreg'}}}
44+
% [loss,model,stats] = bci_train('Data',data, 'Approach','ParadigmFBCSP, 'TargetMarkers',{'n1','n2','n3'})
45+
%
46+
% References;
47+
% [1] Quadrianto Novi, Cuntai Guan, Tran Huy Dat, and Ping Xue, "Sub-band Common Spatial Pattern (SBCSP) for Brain-Computer Interface"
48+
% Proceedings of the 3rd International IEEE EMBS Conference on Neural Engineering Kohala Coast, Hawaii, USA, May 2-5, 2007
49+
% [2] Kai K. Ang, Zhang Y. Chin, Haihong Zhang, Cuntai Guan, "Filter Bank Common Spatial Pattern (FBCSP) in Brain-Computer Interface"
50+
% In 2008 IEEE International Joint Conference on Neural Networks (IEEE World Congress on Computational Intelligence) (June 2008), pp. 2390-2397.
51+
% [3] Owen, A. M., McMillan, K. M., Laird,A. R. & Bullmore, E. "N-back working memory paradigm: A meta-analysis of normative functional neuroimaging studies."
52+
% Human Brain Mapping, 25, 46-59, 2005
53+
%
54+
% Name:
55+
% Filter-Bank CSP
56+
%
57+
% Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
58+
% 2010-04-29
59+
60+
methods
61+
62+
function defaults = preprocessing_defaults(self)
63+
% define the default pre-processing parameters of this paradigm
64+
defaults = {'EpochExtraction',[0.5 3.5],'Resampling',200};
65+
end
66+
67+
function model = feature_adapt(self,varargin)
68+
% adapt a feature representation using the CSP algorithm
69+
args = arg_define(varargin, ...
70+
arg_norep('signal'), ...
71+
arg({'patterns','PatternPairs'},3,uint32([1 1 64 10000]),'CSP patterns per band (times two).','cat','Feature Extraction'), ...
72+
arg({'freqwnds','FreqWindows'},[0.5 3; 4 7; 8 12; 13 30; 31 42],[0 0.5 200 1000],'Frequency bands of interest. Matrix containing one row for the start and end of each frequency band from which CSP patterns shall be computed. Values in Hz.','cat','Feature Extraction'), ...
73+
arg({'winfunc','WindowFunction'},'rect',{'barthann','bartlett','blackman','blackmanharris','bohman','cheb','flattop','gauss','hamming','hann','kaiser','nuttall','parzen','rect','taylor','triang','tukey'},'Type of window function. Typical choices are rect (rectangular), hann, gauss, blackman and kaiser.'),...
74+
arg({'winparam','WindowParameter','param'},[],[],'Parameter of the window function. This is mandatory for cheb, kaiser and tukey and optional for some others.','shape','scalar'),...
75+
arg({'nfft','NFFT'}, [], [],'Size of the FFT used in spectrum calculation. Default value is the greater of 256 or the next power of 2 greater than the length of the signal.' ),...
76+
arg({'winlen','WinLen'},100, [10, 1000], 'Divide the signal into sections of this length for Welch spectrum calculation.'),...
77+
arg({'numoverlap','NumOverlap'}, [], [10, 1000], 'Number of overlap samples from section to next for Welch spectrum calculation.'));
78+
79+
if args.signal.nbchan == 1
80+
error('Multi-band CSP does intrinsically not support single-channel data (it is a spatial filter).'); end
81+
if args.signal.nbchan < args.patterns
82+
error('Multi-band CSP prefers to work on at least as many channels as you request output patterns. Please reduce the number of pattern pairs.'); end
83+
84+
85+
[signal, nof, freqwnds, winfunc, winparam, nfft, winlen, numoverlap] = deal(args.signal, args.patterns, args.freqwnds,...
86+
args.winfunc, args.winparam, args.nfft, args.winlen, args.numoverlap);
87+
88+
[C,S,dum] = size(signal.data);
89+
Fs = signal.srate;
90+
91+
if isempty(numoverlap)
92+
numoverlap = floor(0.5*winlen); end
93+
94+
if (winlen > S) || (numoverlap > winlen)
95+
error(' In Welch method, the length of the window should be smaller than the signal length, and the number of overlap should be smaller than the window length.');
96+
else
97+
win = window_func(winfunc,winlen,winparam);
98+
nwin = floor((S-numoverlap)/(winlen-numoverlap));
99+
end
100+
101+
% The innfft is used internally to design and apply CSP filters
102+
if isempty(nfft)
103+
innfft = 2^(nextpow2(signal.pnts));
104+
else
105+
innfft = max(nfft, 2^(nextpow2(signal.pnts)));
106+
if innfft > nfft
107+
error(' The chosen length of FFT is too short. '); end
108+
end
109+
110+
allfreqs = 0:Fs/innfft:Fs;
111+
allfreqs = allfreqs(1:innfft);
112+
113+
for c=1:2
114+
% compute the per-class epoched data X and its Fourier transform (along time), Xfft
115+
X{c} = exp_eval_optimized(set_picktrials(signal,'rank',c));
116+
[C,S,T] = size(X{c}.data);
117+
Xfft{c} = zeros(C,T,nfft,nwin);
118+
119+
signal_idx = bsxfun(@plus,[1:winlen]',[0:nwin-1]*(winlen-numoverlap));
120+
Xdata_seg = repmat(X{c}.data,1,1,1,nwin); %Xdata_seg -> C,S,T,nwin
121+
Xdata_seg2 = permute(Xdata_seg,[1 3 2 4]); %Xdata_seg2 ->C,T,S,nwin
122+
Xdata_seg3 = reshape(Xdata_seg2(:,:,signal_idx),C,T,[],nwin); % Xdata_seg3 -> C,T,winlen,nwin
123+
Xdata_win = bsxfun(@times,Xdata_seg3,reshape(win,1,1,winlen,1)); % Xdata_win -> C,T,winlen,nwin
124+
Xfft{c} = fft(Xdata_win,innfft,3); % Xfft -> C,T,nfft,nwin
125+
end
126+
127+
filters = [];
128+
patterns = [];
129+
alphas = [];
130+
for fb=1:size(args.freqwnds,1)
131+
132+
[freqs,findx] = getfgrid(Fs,innfft,freqwnds(fb,:));
133+
I = zeros(innfft,1); I(findx)=1;
134+
for cc=1:2
135+
[C,S,T] = size(X{cc}.data);
136+
Xspec{cc} = zeros(C,C);
137+
%Xfft_crop -> C,T,numbands,numwin
138+
Xfft_crop = Xfft{cc}(:, :,findx, :);
139+
F = 2 * real(squeeze(sum(bsxfun(@times,conj(permute(Xfft_crop,[1,5,3,2,4])),permute(Xfft_crop,[5,1,3,2,4])),5))./nwin); % F{c}-> C,C,freqs,T
140+
% compute the cross-spectrum as an average over trials
141+
Xspec{cc} = sum(squeeze(mean(F,4)),3);
142+
end
143+
144+
[V,D] = eig(Xspec{1},Xspec{1}+Xspec{2});
145+
P = inv(V);
146+
filters = [filters V(:,[1:nof end-nof+1:end])];
147+
patterns = [patterns P([1:nof end-nof+1:end],:)'];
148+
alphas = [alphas repmat(I,1, 2*nof)];
149+
150+
end
151+
freq_args.nfft = innfft;
152+
freq_args.win = win;
153+
freq_args.winlen = winlen;
154+
freq_args.numoverlap = numoverlap;
155+
freq_args.nwin = nwin;
156+
model = struct('filters',{filters},'patterns',{patterns},'alphas',{alphas},'freq_args',{freq_args},'chanlocs',{args.signal.chanlocs});
157+
end
158+
159+
function features = feature_extract(self,signal,featuremodel)
160+
freq_args = featuremodel.freq_args;
161+
[nfft, win, winlen, numoverlap, nwin] = deal(freq_args.nfft, freq_args.win, freq_args.winlen, freq_args.numoverlap, freq_args.nwin);
162+
X = signal.data;
163+
[C,S,T] = size(X);
164+
numf = size(featuremodel.filters,2);
165+
features = zeros(T,numf);
166+
167+
Xtemp = permute(X,[3,2,1]); % T,S,C
168+
Xtemp2 = zeros(T,S,size(featuremodel.filters,2));
169+
for t=1:T
170+
Xtemp2(t,:,:) = squeeze(Xtemp(t,:,:)) * featuremodel.filters;
171+
end
172+
%Xtemp2 -> T,S,numf
173+
174+
% Using welch method
175+
signal_idx = bsxfun(@plus,[1:winlen]',[0:nwin-1]*(winlen-numoverlap));
176+
Xdata_seg = repmat(Xtemp2,1,1,1,nwin); %Xdata_seg -> T,S,numf,nwin
177+
Xdata_seg2 = permute(Xdata_seg,[1 3 2 4]); %Xdata_seg2 ->T,numf,S,nwin
178+
Xdata_seg3 = reshape(Xdata_seg2(:,:,signal_idx),T,numf,[],nwin); % Xdata_seg3 -> T,numf,winlen,nwin
179+
Xdata_win = bsxfun(@times,Xdata_seg3,reshape(win,1,1,winlen,1)); % Xdata_win -> T,numf,winlen,nwin
180+
Xfft = fft(Xdata_win,nfft,3); % Xfft -> T,numf,nfft,nwin
181+
Xdata1 = squeeze(sum(Xfft,4)./nwin); % Xdata1 ->T,numf,nfft
182+
Xdata = permute(Xdata1,[1,3,2]); % Xdata ->T,nfft,numf
183+
184+
%Xdata = fft(Xtemp2,nfft,2); %Using regular fft
185+
186+
Xdata2 = bsxfun(@times,Xdata,permute(featuremodel.alphas,[3,1,2]));
187+
Xdata3 = ifft(Xdata2,nfft,2);
188+
Xdata4 = 2*real(Xdata3(:,1:S,:));
189+
features = log(squeeze(var(Xdata4,0,2)));
190+
191+
end
192+
193+
function visualize_model(self,varargin) %#ok<*INUSD>
194+
args = arg_define([0 3],varargin, ...
195+
arg_norep({'myparent','Parent'},[],[],'Parent figure.'), ...
196+
arg_norep({'featuremodel','FeatureModel'},[],[],'Feature model. This is the part of the model that describes the feature extraction.'), ...
197+
arg_norep({'predictivemodel','PredictiveModel'},[],[],'Predictive model. This is the part of the model that describes the predictive mapping.'), ...
198+
arg({'patterns','PlotPatterns'},true,[],'Plot patterns instead of filters. Whether to plot spatial patterns (forward projections) rather than spatial filters.'), ...
199+
arg({'weight_scaled','WeightScaled'},false,[],'Scaled by weight. Whether to scale the patterns by weight.'));
200+
arg_toworkspace(args);
201+
202+
% find the relevant components
203+
scores = predictivemodel.model.w;
204+
scores = sqrt(abs(scores));
205+
% optionally remove the bias if included in w
206+
if length(scores) == size(featuremodel.patterns,2)+1
207+
scores = scores(1:end-1); end
208+
% frequency labels
209+
% titles = repmat({'delta','theta','alpha','beta','gamma'},8,1); titles = titles(:);
210+
% extract relevant patterns
211+
patterns = featuremodel.patterns(:,find(scores)); %#ok<FNDSB>
212+
filters = featuremodel.filters(:,find(scores)); %#ok<FNDSB>
213+
% plot them
214+
if args.weight_scaled
215+
if args.patterns
216+
topoplot_grid(patterns,featuremodel.chanlocs,'scales',scores(find(scores))/max(scores)*1);
217+
else
218+
topoplot_grid(filters,featuremodel.chanlocs,'scales',scores(find(scores))/max(scores)*1);
219+
end
220+
else
221+
if args.patterns
222+
topoplot_grid(patterns,featuremodel.chanlocs);
223+
else
224+
topoplot_grid(filters,featuremodel.chanlocs);
225+
end
226+
end
227+
% figure;
228+
end
229+
230+
function layout = dialog_layout_defaults(self)
231+
% define the default configuration dialog layout
232+
layout = {'SignalProcessing.Resampling.SamplingRate', 'SignalProcessing.EpochExtraction', '', ...
233+
'Prediction.FeatureExtraction.FreqWindows', 'Prediction.FeatureExtraction.TimeWindows', ...
234+
'Prediction.FeatureExtraction.WindowFunction', '', 'Prediction.FeatureExtraction.PatternPairs', '', ...
235+
'Prediction.MachineLearning.Learner'};
236+
end
237+
238+
function tf = needs_voting(self)
239+
tf = true;
240+
end
241+
end
242+
end
243+

0 commit comments

Comments
 (0)