Mercurial > hg > ltpda
comparison m-toolbox/classes/+utils/@math/welchdft.m @ 0:f0afece42f48
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 % WELCHDFT welch method with dft | |
| 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 3 % WELCHDFT comput welch'averaged periodogram with dft | |
| 4 % | |
| 5 % CALL | |
| 6 % Sf = welchdft(x,fs,f,Ns,olap,navs,order,win,psll) | |
| 7 % Sf = welchdft(x,fs,f,Ns,olap,navs,order,win,[]) | |
| 8 % Sf = welchdft(x,fs,f,[],olap,navs,order,win,psll) | |
| 9 % Sf = welchdft(x,fs,f,[],olap,navs,order,win,[]) | |
| 10 % | |
| 11 % | |
| 12 % INPUT | |
| 13 % | |
| 14 % - x, data series, Nx1 double | |
| 15 % - fs, sampling frequency Hz, 1x1 double | |
| 16 % - f, frequenci vector Hz, 1x1 double | |
| 17 % - Ns, length of overlapping segments | |
| 18 % - olap, overlap percentage | |
| 19 % - navs, number of desired averages | |
| 20 % - order, detrend order, -1,0,1,2,3,4,... | |
| 21 % - win, window name. e.g 'BH92' | |
| 22 % - psll, for Kaiser window only | |
| 23 % | |
| 24 % REFERENCES | |
| 25 % | |
| 26 % D. B. Percival and A. T. Walden, Spectral Analysis for Physical | |
| 27 % Applications (Cambridge University Press, Cambridge, 1993) p 291. | |
| 28 % | |
| 29 % L Ferraioli 28-03-2011 | |
| 30 % | |
| 31 % $Id: welchdft.m,v 1.1 2011/03/28 16:37:23 luigi Exp $ | |
| 32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 33 function Sf = welchdft(x,fs,f,Ns,olap,navs,order,win,psll) | |
| 34 | |
| 35 N = numel(x); | |
| 36 if isempty(olap) | |
| 37 % compute shift factor | |
| 38 n = floor((N-Ns)/(navs-1)); | |
| 39 else | |
| 40 if olap>1 | |
| 41 olap = olap/100; | |
| 42 end | |
| 43 % compute Ns | |
| 44 Ns = floor(N/(olap*(navs-1)+1)); | |
| 45 % compute shift factor | |
| 46 n = olap*Ns; | |
| 47 end | |
| 48 | |
| 49 Nf = numel(f); | |
| 50 Sf = zeros(Nf,1); | |
| 51 | |
| 52 % run welch averages | |
| 53 for ii=1:navs | |
| 54 idx1 = n*(ii-1)+1; | |
| 55 idx2 = n*(ii-1)+Ns; | |
| 56 if idx2>N | |
| 57 break | |
| 58 end | |
| 59 seg = x(idx1:idx2); | |
| 60 | |
| 61 % get estimate for a segment | |
| 62 segxx = utils.math.computeDftPeriodogram(seg,fs,f,order,win,psll); | |
| 63 % update sum | |
| 64 Sf = Sf + segxx; | |
| 65 end | |
| 66 % complete average | |
| 67 Sf = Sf./navs; | |
| 68 | |
| 69 | |
| 70 | |
| 71 | |
| 72 | |
| 73 | |
| 74 | |
| 75 | |
| 76 end |