fitlinesc.m 4.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101
  1. function [datafit,Amps,freqs,Fval,sig]=fitlinesc(data,params,p,plt,f0)
  2. % fits significant sine waves to data (continuous data).
  3. %
  4. % Usage: [datafit,Amps,freqs,Fval,sig]=fitlinesc(data,params,p,plt,f0)
  5. %
  6. % Inputs:
  7. % Note that units of Fs, fpass have to be consistent.
  8. % data (data in [N,C] i.e. time x channels/trials or a single
  9. % vector) - required.
  10. % params structure containing parameters - params has the
  11. % following fields: tapers, Fs, fpass, pad
  12. % tapers : precalculated tapers from dpss or in the one of the following
  13. % forms:
  14. % (1) A numeric vector [TW K] where TW is the
  15. % time-bandwidth product and K is the number of
  16. % tapers to be used (less than or equal to
  17. % 2TW-1).
  18. % (2) A numeric vector [W T p] where W is the
  19. % bandwidth, T is the duration of the data and p
  20. % is an integer such that 2TW-p tapers are used. In
  21. % this form there is no default i.e. to specify
  22. % the bandwidth, you have to specify T and p as
  23. % well. Note that the units of W and T have to be
  24. % consistent: if W is in Hz, T must be in seconds
  25. % and vice versa. Note that these units must also
  26. % be consistent with the units of params.Fs: W can
  27. % be in Hz if and only if params.Fs is in Hz.
  28. % The default is to use form 1 with TW=3 and K=5
  29. %
  30. % Fs (sampling frequency) -- optional. Defaults to 1.
  31. % fpass (frequency band to be used in the calculation in the form
  32. % [fmin fmax])- optional.
  33. % Default all frequencies between 0 and Fs/2
  34. % pad (padding factor for the FFT) - optional (can take values -1,0,1,2...).
  35. % -1 corresponds to no padding, 0 corresponds to padding
  36. % to the next highest power of 2 etc.
  37. % e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
  38. % to 512 points, if pad=1, we pad to 1024 points etc.
  39. % Defaults to 0.
  40. % p (P-value to calculate error bars for) - optional.
  41. % Defaults to 0.05/N where N is data length.
  42. % plt (y/n for plot and no plot respectively) - plots the
  43. % Fratio at all frequencies if y
  44. % f0 frequencies at which you want to remove the
  45. % lines - if unspecified the program
  46. % will compute the significant lines
  47. %
  48. %
  49. % Outputs:
  50. % datafit (linear superposition of fitted sine waves)
  51. % Amps (amplitudes at significant frequencies)
  52. % freqs (significant frequencies)
  53. % Fval (Fstatistic at all frequencies)
  54. % sig (significance level for F distribution p value of p)
  55. data=change_row_to_column(data);
  56. [N,C]=size(data);
  57. if nargin < 2 || isempty(params); params=[]; end;
  58. [tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
  59. clear pad fpass err trialave;
  60. if nargin < 3 || isempty(p);p=0.05/N;end;
  61. if nargin < 4 || isempty(plt); plt='n'; end;
  62. if nargin < 5; f0=[]; end;
  63. params.tapers=dpsschk(tapers,N,Fs); % calculate the tapers
  64. [Fval,A,f,sig] = ftestc(data,params,p,plt);
  65. if isempty(f0);
  66. fmax=findpeaks(Fval,sig);
  67. freqs=cell(1,C);
  68. Amps=cell(1,C);
  69. datafit=data;
  70. for ch=1:C;
  71. fsig=f(fmax(ch).loc);
  72. freqs{ch}=fsig;
  73. Amps{ch}=A(fmax(ch).loc,ch);
  74. Nf=length(fsig);
  75. % fprintf('The significant lines for channel %d and the amplitudes are \n',ch);
  76. % for nf=1:Nf;
  77. % fprintf('%12.8f\n',fsig(nf));
  78. % fprintf('%12.8f\n',real(A(fmax(ch).loc(nf),ch)));
  79. % fprintf('%12.8f\n',imag(A(fmax(ch).loc(nf),ch)));
  80. % fprintf('\n');
  81. % end;
  82. datafit(:,ch)=exp(i*2*pi*(0:N-1)'*fsig/Fs)*A(fmax(ch).loc,ch)+exp(-i*2*pi*(0:N-1)'*fsig/Fs)*conj(A(fmax(ch).loc,ch));
  83. end;
  84. else
  85. indx = zeros( length(f0) );
  86. for n=1:length(f0);
  87. [fsig,indx(n)]=min(abs(f-f0(n)));
  88. end;
  89. fsig=f(indx);
  90. for ch=1:C;
  91. freqs{ch}=fsig;
  92. Amps{ch}=A(indx,ch);
  93. Nf=length(fsig);
  94. % fprintf('For channel %d the amplitudes and the Fstatistic at f=%f are \n',ch,f0);
  95. % fprintf('Fstatistic = %12.8f Fthreshold = %12.8f\n',Fval(indx),sig);
  96. % fprintf('Real part of amplitude = %12.8f\n',real(A(indx,ch)));
  97. % fprintf('Imaginary part of amplitude = %12.8f\n',imag(A(indx,ch)));
  98. datafit(:,ch)=exp(i*2*pi*(0:N-1)'*fsig/Fs)*A(indx,ch)+exp(-i*2*pi*(0:N-1)'*fsig/Fs)*conj(A(indx,ch));
  99. end;
  100. end;