Mercurial > hg > ltpda

annotate m-toolbox/classes/@ao/wosa.m @ 49:0bcdf74587d1 database-connection-manager

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Cleanup
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 07 Dec 2011 17:24:36 +0100
parents bc767aaa99a8
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
43
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
1 % WOSA implements Welch's overlaped segmented averaging algorithm with
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
2 % segment detrending and variance estimation.
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
3 %
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
4 % [pxx, f, info] = wosa(x,type,pl)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
5 % [pxx, f, info] = wosa(x,y,type,pl)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
6 %
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
7 % INPUTS: x - input analysis objects
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
8 % y - input analysis objects
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
9 % type - type of estimation:
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
10 % 'psd' - compute Power Spectral Denstiy (PSD)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
11 % 'cpsd' - compute cross-spectral density
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
12 % 'tfe' - estimate transfer function between inputs
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
13 % 'mscohere' - estimate magnitude-squared cross-coherence
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
14 % 'cohere' - estimate complex cross-coherence
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
15 % pl - input parameter list
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
16 %
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
17 % PARAMETERS: 'Win' - a specwin window object [default: Kaiser -200dB psll]
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
18 % 'Olap' - segment percent overlap [default: taken from window function]
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
19 % 'Nfft' - number of samples in each fft [default: length of input data]
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
20 % 'Scale' - one of
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
21 % 'ASD' - amplitude spectral density
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
22 % 'PSD' - power spectral density [default]
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
23 % 'AS' - amplitude spectrum
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
24 % 'PS' - power spectrum
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
25 % * applies only to spectrum 'Type' 'psd'
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
26 % 'Order' - order of segment detrending
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
27 % -1 - no detrending
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
28 % 0 - subtract mean [default]
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
29 % 1 - subtract linear fit
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
30 % N - subtract fit of polynomial, order N
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
31 %
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
32 % Version: $Id: wosa.m,v 1.5 2011/12/02 07:08:11 hewitson Exp $
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
33 %
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
34
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
35 function varargout = wosa(varargin)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
36 import utils.const.*
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
37
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
38 % Process inputs
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
39 if nargin == 3
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
40 a = varargin{1};
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
41 esttype = varargin{2};
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
42 pl = varargin{3};
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
43 inunits = a.data.yunits;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
44 L = a.len;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
45 else
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
46 a = varargin{1};
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
47 b = varargin{2};
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
48 esttype = varargin{3};
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
49 pl = varargin{4};
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
50 if a.data.fs ~= b.data.fs
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
51 error('The two time-series have different sample rates.');
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
52 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
53 inunits = b.data.yunits / a.data.yunits;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
54 L = min(a.len, b.len);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
55 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
56
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
57 % Parse inputs
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
58 win = find(pl, 'Win');
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
59 nfft = find(pl, 'Nfft');
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
60 percentOlap = find(pl, 'Olap')/100;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
61 scale = find(pl, 'scale');
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
62 xOlap = round(percentOlap*nfft);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
63 detrendOrder = find(pl, 'order');
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
64 fs = a.fs;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
65 winVals = win.win.'; % because we always get a column from ao.y
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
66
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
67 % Compute segment details
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
68
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
69 nSegments = fix((L - xOlap)./(nfft - xOlap));
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
70 utils.helper.msg(msg.PROC3, 'N segment: %d', nfft);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
71 utils.helper.msg(msg.PROC3, 'N overlap: %d', xOlap);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
72 utils.helper.msg(msg.PROC3, 'N segments: %d', nSegments);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
73
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
74 % Compute start and end indices of each segment
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
75 segmentStep = nfft-xOlap;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
76 segmentStarts = 1:segmentStep:nSegments*segmentStep;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
77 segmentEnds = segmentStarts+nfft-1;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
78
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
79 % Estimate the averaged periodogram for the desired quantity
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
80 switch lower(esttype)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
81 case 'psd'
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
82 % Compute averaged periodogram
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
83 [Sxx, Svxx] = psdPeriodogram(a, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
84 case 'cpsd'
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
85 [Sxx, Svxx] = cpsdPeriodogram(a, b, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
86 case 'tfe'
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
87 [Sxx, Sxy, Syy] = tfePeriodogram(a, b, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
88 case {'mscohere','cohere'}
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
89 [Sxx, Sxy, Syy] = tfePeriodogram(a, b, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
90 otherwise
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
91 error('Unknown estimation type %s', esttype);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
92 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
93
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
94 % Scale to PSD
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
95 switch lower(esttype)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
96 case {'psd','cpsd'}
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
97 [P, Pvxx] = scaleToPSD(Sxx, Svxx, nfft, fs);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
98 % the 1/nSegments factor should come after welchscale if we don't
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
99 % want to apply sqrt() to it.
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
100 % We correct for that here. It is only needed for 'asd','as' in
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
101 % psd/cpsd, the other cases go always through 'PSD'.
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
102 if (strcmpi(scale,'PSD') || strcmpi(scale,'PS'))
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
103 dP = Pvxx;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
104 elseif (strcmpi(scale,'ASD') || strcmpi(scale,'AS'))
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
105 dP = Pvxx/nSegments;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
106 else
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
107 error('### Unknown scale')
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
108 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
109 case 'tfe'
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
110 % Compute the 1-sided or 2-sided PSD [Power/freq] or mean-square [Power].
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
111 % Also, corresponding freq vector and freq units.
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
112 % In the Cross PSD, the frequency vector and xunits are not used.
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
113 Pxx = scaleToPSD(Sxx, [], nfft, fs);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
114 Pxy = scaleToPSD(Sxy, [], nfft, fs);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
115 Pyy = scaleToPSD(Syy, [], nfft, fs);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
116 % mean and std
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
117 P = Pxy ./ Pxx; % Txy
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
118 if nSegments == 1
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
119 dP =[];
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
120 else
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
121 dP = (nSegments/(nSegments-1)^2)*(Pyy./Pxx).*(1 - (abs(Pxy).^2)./(Pxx.*Pyy));
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
122 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
123 case 'mscohere'
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
124 % Magnitude Square Coherence estimate.
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
125 % Auto PSD for 2nd input vector. The freq vector & xunits are not
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
126 % used.
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
127 Pxx = scaleToPSD(Sxx, [], nfft, fs);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
128 Pxy = scaleToPSD(Sxy, [], nfft, fs);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
129 Pyy = scaleToPSD(Syy, [], nfft, fs);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
130 % mean and std
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
131 P = (abs(Pxy).^2)./(Pxx.*Pyy); % Magnitude-squared coherence
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
132 dP = (2*P/nSegments).*(1-P).^2;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
133 case 'cohere'
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
134 % Complex Coherence estimate.
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
135 % Auto PSD for 2nd input vector. The freq vector & xunits are not
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
136 % used.
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
137 Pxx = scaleToPSD(Sxx, [], nfft, fs);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
138 Pxy = scaleToPSD(Sxy, [], nfft, fs);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
139 Pyy = scaleToPSD(Syy, [], nfft, fs);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
140 P = Pxy./sqrt(Pxx.*Pyy); % Complex coherence
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
141 dP = (2*abs(P)/nSegments).*(1-abs(P)).^2;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
142
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
143 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
144
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
145 % Compute frequencies
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
146 freqs = psdfreqvec('npts', nfft,'Fs', fs, 'Range', 'half').';
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
147
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
148 % Scale to required units
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
149 [Pxx, dP, info] = utils.math.welchscale(P, dP, winVals, fs, scale, inunits);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
150 info.navs = nSegments;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
151
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
152 if nSegments ==1
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
153 dev = [];
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
154 else
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
155 dev = sqrt(dP);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
156 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
157
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
158 % Set outputs
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
159 varargout = {Pxx, freqs, info, dev};
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
160
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
161 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
162
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
163 % scale averaged periodogram to PSD
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
164 function [Pxx, Pvxx] = scaleToPSD(Sxx, Svxx, nfft, fs)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
165
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
166 % Take 1-sided spectrum which means we double the power in the
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
167 % appropriate bins
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
168 if rem(nfft,2),
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
169 indices = 1:(nfft+1)/2; % ODD
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
170 Sxx1sided = Sxx(indices,:);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
171 % double the power except for the DC bin
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
172 Sxx = [Sxx1sided(1,:); 2*Sxx1sided(2:end,:)];
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
173 if ~isempty(Svxx)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
174 Svxx1sided = Svxx(indices,:);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
175 Svxx = [Svxx1sided(1,:); 4*Svxx1sided(2:end,:)];
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
176 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
177 else
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
178 indices = 1:nfft/2+1; % EVEN
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
179 Sxx1sided = Sxx(indices,:);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
180 % Double power except the DC bin and the Nyquist bin
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
181 Sxx = [Sxx1sided(1,:); 2*Sxx1sided(2:end-1,:); Sxx1sided(end,:)];
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
182 if ~isempty(Svxx)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
183 Svxx1sided = Svxx(indices,:); % Take only [0,pi] or [0,pi)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
184 Svxx = [Svxx1sided(1,:); 4*Svxx1sided(2:end-1,:); Svxx1sided(end,:)];
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
185 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
186 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
187
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
188 % Now scale to PSD
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
189 Pxx = Sxx./fs;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
190 Pvxx = Svxx./fs^2;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
191
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
192 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
193
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
194 % compute tfe
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
195 function [Sxx, Sxy, Syy] = tfePeriodogram(x, y, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
196
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
197 nfft = segmentEnds(1);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
198 Sxx = zeros(nfft,1); % Initialize Sxx
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
199 Sxy = zeros(nfft,1); % Initialize Sxy
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
200 Syy = zeros(nfft,1); % Initialize Syy
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
201 % loop over segments
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
202 for ii = 1:nSegments
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
203 if detrendOrder < 0
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
204 xseg = x.y(segmentStarts(ii):segmentEnds(ii));
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
205 yseg = y.y(segmentStarts(ii):segmentEnds(ii));
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
206 else
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
207 [xseg,coeffs] = ltpda_polyreg(x.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
208 [yseg,coeffs] = ltpda_polyreg(y.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
209 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
210
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
211 % Compute periodograms
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
212 Sxxk = wosa_periodogram(xseg, [], winVals, nfft);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
213 Sxyk = wosa_periodogram(yseg, xseg, winVals, nfft);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
214 Syyk = wosa_periodogram(yseg, [], winVals, nfft);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
215
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
216 Sxx = Sxx + Sxxk;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
217 Sxy = Sxy + Sxyk;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
218 Syy = Syy + Syyk;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
219 % don't need to be divided by k because only rations are used here
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
220 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
221
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
222 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
223
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
224 % compute cpsd
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
225 function [Sxx, Svxx] = cpsdPeriodogram(x, y, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
226
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
227 Mnxx = 0;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
228 Mn2xx = 0;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
229 nfft = segmentEnds(1);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
230 for ii = 1:nSegments
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
231 if detrendOrder < 0
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
232 xseg = x.y(segmentStarts(ii):segmentEnds(ii));
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
233 yseg = y.y(segmentStarts(ii):segmentEnds(ii));
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
234 else
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
235 [xseg,coeffs] = ltpda_polyreg(x.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
236 [yseg,coeffs] = ltpda_polyreg(y.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
237 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
238
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
239 % Compute periodogram
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
240 Sxxk = wosa_periodogram(xseg, yseg, winVals, nfft);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
241
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
242 % Welford's algorithm to update mean and variance
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
243 Qxx = Sxxk - Mnxx;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
244 Mnxx = Mnxx +Qxx/ii;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
245 Mn2xx = Mn2xx + abs(Qxx.*conj(Sxxk - Mnxx));
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
246 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
247 Sxx = Mnxx;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
248 if nSegments ==1
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
249 Svxx = [];
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
250 else
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
251 Svxx = Mn2xx/(nSegments-1)/nSegments;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
252 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
253
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
254
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
255 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
256
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
257 % compute psd
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
258 function [Sxx, Svxx] = psdPeriodogram(x, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
259
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
260 Mnxx = 0;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
261 Mn2xx = 0;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
262 nfft = segmentEnds(1);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
263 % Loop over the segments
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
264 for ii = 1:nSegments
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
265 % Detrend if desired
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
266 if detrendOrder < 0
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
267 seg = x.y(segmentStarts(ii):segmentEnds(ii));
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
268 else
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
269 [seg,coeffs] = ltpda_polyreg(x.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
270 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
271 % Compute periodogram
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
272 Sxxk = wosa_periodogram(seg, [], winVals,nfft);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
273 % Welford's algorithm for updating mean and variance
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
274 if ii == 1
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
275 Mnxx = Sxxk;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
276 else
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
277 Qxx = Sxxk - Mnxx;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
278 Mnxx = Mnxx + Qxx/ii;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
279 Mn2xx = Mn2xx + Qxx.*(Sxxk - Mnxx);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
280 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
281 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
282 Sxx = Mnxx;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
283 if nSegments == 1
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
284 Svxx = [];
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
285 else
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
286 Svxx = Mn2xx/(nSegments-1)/nSegments;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
287 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
288
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
289 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
290
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
291 % Scaled periodogram of one or two input signals
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
292 function Sxx = wosa_periodogram(x, y, win, nfft)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
293
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
294 % window data
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
295 xwin = x.*win;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
296 isCross = false;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
297 if ~isempty(y)
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
298 ywin = y.*win;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
299 isCross = true;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
300 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
301
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
302 % take fft
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
303 X = fft(xwin, nfft);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
304 if isCross
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
305 Y = fft(ywin, nfft);
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
306 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
307
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
308 % Compute scale factor to compensate for the window power
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
309 K = win'*win;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
310
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
311 % Compute scaled power
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
312 Sxx = X.*conj(X)/K;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
313 if isCross,
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
314 Sxx = X.*conj(Y)/K;
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
315 end
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
316
bc767aaa99a8 CVS Update
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff changeset
317 end