fits significant sine waves to data (continuous data) using overlapping windows.
Usage: [datac,datafit]=rmlinesmovingwinc(data,movingwin,tau,params,p,plt)
Inputs:
Note that units of Fs, fpass have to be consistent.
data (data in [N,C] i.e. time x channels/trials or as a single vector) - required.
movingwin (in the form [window winstep] i.e length of moving
window and step size)
Note that units here have
to be consistent with
units of Fs - required
tau parameter controlling degree of smoothing for the amplitudes - we use the
function 1-1/(1+exp(-tau*(x-Noverlap/2)/Noverlap) in the region of overlap to smooth
the sinewaves across the overlap region. Noverlap is the number of points
in the overlap region. Increasing tau leads to greater overlap smoothing,
typically specifying tau~10 or higher is reasonable. tau=1 gives an almost
linear smoothing function. tau=100 gives a very steep sigmoidal. The default is tau=10.
params structure containing parameters - params has the
following fields: tapers, Fs, fpass, pad
tapers : precalculated tapers from dpss or in the one of the following
forms:
(1) A numeric vector [TW K] where TW is the
time-bandwidth product and K is the number of
tapers to be used (less than or equal to
2TW-1).
(2) A numeric vector [W T p] where W is the
bandwidth, T is the duration of the data and p
is an integer such that 2TW-p tapers are used. In
this form there is no default i.e. to specify
the bandwidth, you have to specify T and p as
well. Note that the units of W and T have to be
consistent: if W is in Hz, T must be in seconds
and vice versa. Note that these units must also
be consistent with the units of params.Fs: W can
be in Hz if and only if params.Fs is in Hz.
The default is to use form 1 with TW=3 and K=5
Note that T has to be equal to movingwin(1).
Fs (sampling frequency) -- optional. Defaults to 1.
fpass (frequency band to be used in the calculation in the form
[fmin fmax])- optional.
Default all frequencies between 0 and Fs/2
pad (padding factor for the FFT) - optional (can take values -1,0,1,2...).
-1 corresponds to no padding, 0 corresponds to padding
to the next highest power of 2 etc.
e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
to 512 points, if pad=1, we pad to 1024 points etc.
Defaults to 0.
p (P-value to calculate error bars for) - optional.
Defaults to 0.05/Nwin where Nwin is length of window which
corresponds to a false detect probability of approximately 0.05.
plt (y/n for plot and no plot respectively) - default no
plot.
f0 frequencies at which you want to remove the
lines - if unspecified the program uses the f statistic
to determine appropriate lines.
Outputs:
datafit (fitted sine waves)
datac (cleaned up data)0001 function [datac,datafit,Amps,freqs]=rmlinesmovingwinc(data,movingwin,tau,params,p,plt,f0) 0002 % fits significant sine waves to data (continuous data) using overlapping windows. 0003 % 0004 % Usage: [datac,datafit]=rmlinesmovingwinc(data,movingwin,tau,params,p,plt) 0005 % 0006 % Inputs: 0007 % Note that units of Fs, fpass have to be consistent. 0008 % data (data in [N,C] i.e. time x channels/trials or as a single vector) - required. 0009 % movingwin (in the form [window winstep] i.e length of moving 0010 % window and step size) 0011 % Note that units here have 0012 % to be consistent with 0013 % units of Fs - required 0014 % tau parameter controlling degree of smoothing for the amplitudes - we use the 0015 % function 1-1/(1+exp(-tau*(x-Noverlap/2)/Noverlap) in the region of overlap to smooth 0016 % the sinewaves across the overlap region. Noverlap is the number of points 0017 % in the overlap region. Increasing tau leads to greater overlap smoothing, 0018 % typically specifying tau~10 or higher is reasonable. tau=1 gives an almost 0019 % linear smoothing function. tau=100 gives a very steep sigmoidal. The default is tau=10. 0020 % params structure containing parameters - params has the 0021 % following fields: tapers, Fs, fpass, pad 0022 % tapers : precalculated tapers from dpss or in the one of the following 0023 % forms: 0024 % (1) A numeric vector [TW K] where TW is the 0025 % time-bandwidth product and K is the number of 0026 % tapers to be used (less than or equal to 0027 % 2TW-1). 0028 % (2) A numeric vector [W T p] where W is the 0029 % bandwidth, T is the duration of the data and p 0030 % is an integer such that 2TW-p tapers are used. In 0031 % this form there is no default i.e. to specify 0032 % the bandwidth, you have to specify T and p as 0033 % well. Note that the units of W and T have to be 0034 % consistent: if W is in Hz, T must be in seconds 0035 % and vice versa. Note that these units must also 0036 % be consistent with the units of params.Fs: W can 0037 % be in Hz if and only if params.Fs is in Hz. 0038 % The default is to use form 1 with TW=3 and K=5 0039 % Note that T has to be equal to movingwin(1). 0040 % 0041 % Fs (sampling frequency) -- optional. Defaults to 1. 0042 % fpass (frequency band to be used in the calculation in the form 0043 % [fmin fmax])- optional. 0044 % Default all frequencies between 0 and Fs/2 0045 % pad (padding factor for the FFT) - optional (can take values -1,0,1,2...). 0046 % -1 corresponds to no padding, 0 corresponds to padding 0047 % to the next highest power of 2 etc. 0048 % e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT 0049 % to 512 points, if pad=1, we pad to 1024 points etc. 0050 % Defaults to 0. 0051 % p (P-value to calculate error bars for) - optional. 0052 % Defaults to 0.05/Nwin where Nwin is length of window which 0053 % corresponds to a false detect probability of approximately 0.05. 0054 % plt (y/n for plot and no plot respectively) - default no 0055 % plot. 0056 % f0 frequencies at which you want to remove the 0057 % lines - if unspecified the program uses the f statistic 0058 % to determine appropriate lines. 0059 % 0060 % Outputs: 0061 % datafit (fitted sine waves) 0062 % datac (cleaned up data) 0063 if nargin < 2; error('Need data and window parameters'); end; 0064 if nargin < 4 || isempty(params); params=[]; end; 0065 0066 if length(params.tapers)==3 & movingwin(1)~=params.tapers(2); 0067 error('Duration of data in params.tapers is inconsistent with movingwin(1), modify params.tapers(2) to proceed') 0068 end 0069 0070 [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params); % set defaults for params 0071 clear err trialave 0072 if nargin < 6; plt='n'; end; 0073 % 0074 % Window,overlap and frequency information 0075 % 0076 data=change_row_to_column(data); 0077 [N,C]=size(data); 0078 Nwin=round(Fs*movingwin(1)); % number of samples in window 0079 Nstep=round(movingwin(2)*Fs); % number of samples to step through 0080 Noverlap=Nwin-Nstep; % number of points in overlap 0081 % 0082 % Sigmoidal smoothing function 0083 % 0084 if nargin < 3 || isempty(tau); tau=10; end; % smoothing parameter for sigmoidal overlap function 0085 x=(1:Noverlap)'; 0086 smooth=1./(1+exp(-tau.*(x-Noverlap/2)/Noverlap)); % sigmoidal function 0087 smooth=repmat(smooth,[1 C]); 0088 % 0089 % Start the loop 0090 % 0091 if nargin < 5 || isempty(p); p=0.05/Nwin; end % default for p value 0092 if nargin < 7 || isempty(f0); f0=[]; end; % empty set default for f0 - uses F statistics to determine the frequencies 0093 params.tapers=dpsschk(tapers,Nwin,Fs); % check tapers 0094 winstart=1:Nstep:N-Nwin+1; 0095 nw=length(winstart); 0096 datafit=zeros(winstart(nw)+Nwin-1,C); 0097 Amps=cell(1,nw); 0098 freqs=cell(1,nw); 0099 for n=1:nw; 0100 indx=winstart(n):winstart(n)+Nwin-1; 0101 datawin=data(indx,:); 0102 [datafitwin,as,fs]=fitlinesc(datawin,params,p,'n',f0); 0103 Amps{n}=as; 0104 freqs{n}=fs; 0105 datafitwin0=datafitwin; 0106 if n>1; datafitwin(1:Noverlap,:)=smooth.*datafitwin(1:Noverlap,:)+(1-smooth).*datafitwin0(Nwin-Noverlap+1:Nwin,:);end; 0107 datafit(indx,:)=datafitwin; 0108 end; 0109 datac=data(1:size(datafit,1),:)-datafit; 0110 if strcmp(plt,'y'); 0111 [S,f]=mtspectrumsegc(data,movingwin(1),params); 0112 [Sc,fc]=mtspectrumsegc(datac,movingwin(1),params); 0113 plot(f,10*log10(S),fc,10*log10(Sc)); 0114 end;