In [1]:
warning('off','all');
if ~exist('run_once','var')
curdir = pwd;
cd('../..');
startup;
cd(curdir);
run_once = true;
end
In [2]:
%plot -s 600,400
In [3]:
nchannels = 2;
ntrials = 10;
norder = 3;
fs = 250;
nsamples_signal = fs/2;
% generate biosignal
biosignal = BioSignal.muscle('fs',fs,'nsamples',nsamples_signal,'bw',5,'fo',5);
% create erp
erp = [-1*fliplr(biosignal(2:end)) biosignal];
% normalize erp
erp = erp./max(erp);
figure;
nsamples_signal = length(biosignal);
time_start = 0;
time_end = (nsamples_signal*2-1)/fs;
time_axis = linspace(time_start,time_end,nsamples_signal*2-1);
plot(time_axis, erp);
xlim([time_axis(1), time_axis(end)]);
title('ERP');
In [4]:
nsamples = length(erp);
% allocate mem for data
x_orig = zeros(nchannels,nsamples,ntrials);
x = zeros(nchannels,nsamples,ntrials);
for i=1:ntrials
% simulate noise
sigma = 0.05;
noise = mvnrnd([0 0], [sigma 0; 0 sigma],nsamples)';
% add erp variability
H = [2; -5];
[signal_trial,scale] = BioSignal.rand_scale(erp,1,0.05);
%fprintf('scale: %f\n',scale);
[signal_trial,shift] = BioSignal.rand_time_shift(signal_trial,0,15);
%fprintf('shift: %d\n',shift);
% create the signal
x_orig(:,:,i) = H*signal_trial + noise;
% normalize variance of each channel to unit variance
x(:,:,i) = x_orig(:,:,i)./repmat(std(x_orig(:,:,i),0,2),1,nsamples);
end
In [5]:
figure;
plot(time_axis,squeeze(x_orig(1,:,:)));
xlim([time_axis(1), time_axis(end)]);
title('Trials');
In [6]:
figure;
plot(time_axis,squeeze(x(1,:,:)));
xlim([time_axis(1), time_axis(end)]);
title('Trials Normalized');
In [7]:
% take the average of all trials
x_avg = mean(x,3);
figure;
plot(time_axis,x_avg(1,:));
xlim([time_axis(1), time_axis(end)]);
title('Trials Averaged');
In [8]:
plot_options = {'ch1',2,'ch2',2};
In [9]:
order_est = 3;
verbosity = 0;
lambda = 0.99;
filter = MCMTQRDLSL1(ntrials,nchannels,order_est,lambda);
trace = LatticeTrace(filter,'fields',{'Kf'});
% run the filter
figure;
plot_options_cust = [plot_options {'title','Multi Trial - 10 Trials'}];
trace.run(x,'verbosity',verbosity,'mode','plot',...
'plot_options',plot_options_cust);
In [10]:
lambda = 0.99;
filter = MCMTQRDLSL1(2,nchannels,order_est,lambda);
trace = LatticeTrace(filter,'fields',{'Kf'});
% run the filter
figure;
plot_options_cust = [plot_options {'title','Multi Trial - 2 Trials'}];
trace.run(x(:,:,1:2),'verbosity',verbosity,'mode','plot',...
'plot_options',plot_options_cust);
In [11]:
filter = MQRDLSL1(nchannels,order_est,lambda);
% filter = MQRDLSL2(nchannels,order_est,lambda);
trace = LatticeTrace(filter,'fields',{'Kf'});
% run the filter
figure;
plot_options_cust = [plot_options {'title',sprintf('Single Trial - Trial %d',i)}];
trace.run(x(:,:,1),'verbosity',verbosity,'mode','plot',...
'plot_options',plot_options_cust);
In [12]:
filter = MQRDLSL1(nchannels,order_est,lambda);
% filter = MQRDLSL2(nchannels,order_est,lambda);
trace = LatticeTrace(filter,'fields',{'Kf'});
% run the filter
figure;
plot_options_cust = [plot_options {'title','Single Trial - Averaged ERP'}];
trace.run(x_avg,'verbosity',verbosity,'mode','plot',...
'plot_options',plot_options_cust);