+
−% COMPUTEDFTPERIODOGRAM compute periodogram with dft
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−% 
+
−% DESCRIPTION
+
−% 
+
−% Define one sided periodogram as:
+
−% Xp(f) = 2*T*|DFT(w(t)*x(t))/T|^2
+
−% Where w(t) are window sameples, which are supposed to be square
+
−% integrable: sum(w(t)^2)=1
+
−% 
+
−% CALL 
+
−% Sf = computeDftPeriodogram(x,fs,f,order,win,psll)
+
−% Sf = computeDftPeriodogram(x,fs,f,order,win,[])
+
−% 
+
−% 
+
−% INPUT
+
−% 
+
−% - x, data series, Nx1 double
+
−% - fs, sampling frequency Hz, 1x1 double
+
−% - f, frequenci vector Hz, 1x1 double
+
−% - order, detrend order, -1,0,1,2,3,4,...
+
−% - win, window name. e.g 'BH92'
+
−% - psll, for Kaiser window only
+
−% 
+
−% REFERENCES
+
−% 
+
−% D. B. Percival and A. T. Walden, Spectral Analysis for Physical
+
−% Applications (Cambridge University Press, Cambridge, 1993) p 291.
+
−% 
+
−% L Ferraioli 28-03-2011
+
−%
+
−% $Id: computeDftPeriodogram.m,v 1.2 2011/03/29 15:36:43 luigi Exp $
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−function Sf = computeDftPeriodogram(x,fs,f,order,win,psll)
+
−
+
−  
+
−  % detrend segment
+
−  if order < 0
+
−    xd = x; % no detrend
+
−  else
+
−    [xd,coeffs] = ltpda_polyreg(x,order);
+
−  end
+
−  
+
−  N = numel(x);
+
−  
+
−  % get window
+
−  switch lower(win)
+
−    case 'kaiser'
+
−      wnd = specwin('Kaiser', N, psll);
+
−    otherwise
+
−      wnd = specwin(win, N);
+
−  end
+
−  w = wnd.win;
+
−  w = w./sqrt(wnd.ws2);  % compensates for the power of the window.
+
−  
+
−  if size(w)~=size(xd)
+
−    w = w.';
+
−  end
+
−  
+
−  % apply window
+
−  xdw = xd.*w;
+
−  
+
−  T = 1/fs;
+
−  
+
−  Nf = numel(f);
+
−  Sf = zeros(Nf,1);
+
−  for ii=1:Nf
+
−    xx = utils.math.dft(xdw,f(ii),T);
+
−    xx = xx/T;
+
−    Sf(ii) = 2*T*xx*conj(xx);
+
−  end
+
−
+
−  
+
−
+
−end