view m-toolbox/classes/+utils/@math/computeDftPeriodogram.m @ 21:8be9deffe989
database-connection-manager
Update ltpda_uo.update
author
Daniele Nicolodi <nicolodi@science.unitn.it>
date
Mon, 05 Dec 2011 16:20:06 +0100 (2011-12-05)
parents
f0afece42f48
children
line source
+ − % 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