Mercurial > hg > ltpda

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % COMPUTEDFTPERIODOGRAM compute periodogram with dft
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION
5 %
6 % Define one sided periodogram as:
7 % Xp(f) = 2*T*|DFT(w(t)*x(t))/T|^2
8 % Where w(t) are window sameples, which are supposed to be square
9 % integrable: sum(w(t)^2)=1
10 %
11 % CALL
12 % Sf = computeDftPeriodogram(x,fs,f,order,win,psll)
13 % Sf = computeDftPeriodogram(x,fs,f,order,win,[])
14 %
15 %
16 % INPUT
17 %
18 % - x, data series, Nx1 double
19 % - fs, sampling frequency Hz, 1x1 double
20 % - f, frequenci vector Hz, 1x1 double
21 % - order, detrend order, -1,0,1,2,3,4,...
22 % - win, window name. e.g 'BH92'
23 % - psll, for Kaiser window only
24 %
25 % REFERENCES
26 %
27 % D. B. Percival and A. T. Walden, Spectral Analysis for Physical
28 % Applications (Cambridge University Press, Cambridge, 1993) p 291.
29 %
30 % L Ferraioli 28-03-2011
31 %
32 % $Id: computeDftPeriodogram.m,v 1.2 2011/03/29 15:36:43 luigi Exp $
33 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34 function Sf = computeDftPeriodogram(x,fs,f,order,win,psll)
35
36
37 % detrend segment
38 if order < 0
39 xd = x; % no detrend
40 else
41 [xd,coeffs] = ltpda_polyreg(x,order);
42 end
43
44 N = numel(x);
45
46 % get window
47 switch lower(win)
48 case 'kaiser'
49 wnd = specwin('Kaiser', N, psll);
50 otherwise
51 wnd = specwin(win, N);
52 end
53 w = wnd.win;
54 w = w./sqrt(wnd.ws2); % compensates for the power of the window.
55
56 if size(w)~=size(xd)
57 w = w.';
58 end
59
60 % apply window
61 xdw = xd.*w;
62
63 T = 1/fs;
64
65 Nf = numel(f);
66 Sf = zeros(Nf,1);
67 for ii=1:Nf
68 xx = utils.math.dft(xdw,f(ii),T);
69 xx = xx/T;
70 Sf(ii) = 2*T*xx*conj(xx);
71 end
72
73
74
75 end