vb_polysincosfilt

 A kind of high pass filter

[Y,A,Base] = vb_polysincosfilt(X,basis,sampHz,flag_const)

Y : drift free data

X : Ntrial(Nsensor) * Nsamp
basis : sine wave cycle  e.g. [4 8 16]
sampHz:  e.g. 1000
flag_const: 1 - add constant term in basis matrix.


 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [Y,A,Base] = vb_polysincosfilt(X,basis,sampHz,flag_const)
0002 % A kind of high pass filter
0003 %
0004 %[Y,A,Base] = vb_polysincosfilt(X,basis,sampHz,flag_const)
0005 %
0006 %Y : drift free data
0007 %
0008 %X : Ntrial(Nsensor) * Nsamp
0009 %basis : sine wave cycle  e.g. [4 8 16]
0010 %sampHz:  e.g. 1000
0011 %flag_const: 1 - add constant term in basis matrix.
0012 %
0013 %
0014 % Copyright (C) 2011, ATR All Rights Reserved.
0015 % License : New BSD License(see VBMEG_LICENSE.txt)
0016 
0017 if nargin<4
0018   flag_const = 1;
0019 end
0020 
0021 isi = 1/sampHz;
0022 
0023 [Nsensor,Nsamp] = size(X);
0024 t = isi:isi:isi*Nsamp;
0025 
0026 Nbase = length(basis);
0027 
0028 for i = 1:Nbase
0029     Base(i,:) = sin(2*pi*t*(1/basis(i)));
0030     Base(Nbase+i,:) = cos(2*pi*t*(1/basis(i)));
0031 end
0032 if flag_const
0033   Base(Nbase*2+1,:) = 1;
0034 end
0035 
0036 %model : X = A*Base
0037 A = X*Base'*pinv(Base*Base');
0038 Drift = A*Base;
0039 Y = X - Drift;