diff m-toolbox/classes/+utils/@math/computeDftPeriodogram.m @ 0:f0afece42f48

Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/computeDftPeriodogram.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,75 @@
+% 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
\ No newline at end of file