Mercurial > hg > ltpda

comparison m-toolbox/m/sigproc/frequency_domain/ltpda_spsd.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 % ltpda_spsd smooths a spectrum.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: ltpda_spsd smooths a spectrum.
5 %
6 % CALL: [freqsOut, pow, nFreqs, sigmaP, sumMat] = ltpda_spsd(freqs, pow, linCoef, logCoef, sumMat,nFreqs)
7 %
8 % INPUTS: freqsOut - frequency vector, can be left empty
9 % pow - power spectrum (density)
10 % nFreqs - frequency intervals used to derive averages
11 % linCoef, logCoef
12 % - values to use to scale the smoothing
13 % averaging filter. It will be in linCoef.NBins^logCoef.
14 % logCoef should be 2/3 for PSD and 4/5 for PSD data
15 % which is later plotted in ASD scale.
16 % mode - 'keepFrequencies' convoles using a filter, leaving all
17 % frequencies but they are correlated
18 % 'removeFrequencies' convoles using a filter, leaving
19 % only data at some uncorrelated frequencies
20 % 'addUp' like above, but sums up instead of takin a
21 % mean value.
22 %
23 % OUTPUTS: freqs, pow : output frequency and power vector.
24 %
25 % <a href="matlab:web(ao.getInfo('spsd').tohtml, '-helpbrowser')">Parameter Sets</a>
26 %
27 % VERSION: $Id: ltpda_spsd.m,v 1.4 2011/02/21 16:55:21 adrien Exp $
28 %
29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 function varargout = ltpda_spsd(varargin)
31 %% collecting inputs
32 freqs = varargin{1};
33 pow = varargin{2};
34 L = length(pow);
35 pow = reshape(pow, [L, 1]);
36 freqs = reshape(freqs, [numel(freqs), 1]);
37 linCoef = varargin{3};
38 logCoef = varargin{4};
39 if nargin >4
40 % optional inputs include the averaging matrix
41 sumMat = varargin{5};
42 nFreqs = varargin{6};
43 else
44 % averaging matrix must be computed from input data
45 [nFreqs, sumMat] = getMatrix(L, linCoef, logCoef);
46 end
47 %% evaluating output
48 computeVar = nargout>=4;
49 computeFreqs = ~isempty(freqs);
50 sumPow = sumMat * pow; % power sums
51 powOut = sumPow ./ nFreqs; % power average
52 if computeFreqs
53 freqsOut = (sumMat * freqs) ./ nFreqs; % correponding mean frequency
54 else
55 freqsOut = [];
56 end
57 if computeVar
58 nBinsEff = sumPow.^2 ./ (sumMat * pow.^2) ; % number of Chi2 variables for STD
59 else
60 nBinsEff = [];
61 end
62
63 %% allocating output
64 varargout = {freqsOut, powOut, nFreqs, nBinsEff, sumMat};
65
66 end
67
68
69 function [nFreqs, sumMat] = getMatrix(L, linCoef, logCoef)
70 %% initializing loop
71 iAvg = 1;
72 iMin = 1;
73 idxAvg = zeros(1,L); % index of corresponding average
74 widths = ceil(linCoef):ceil(linCoef * L^logCoef + 2); % proposed width of averagin filter
75 iMaxForWidths = min(L, floor( (widths./linCoef).^(1/logCoef) + widths-1 ) ); % last sample to be processed with a given filter width
76 lastWidth = find(iMaxForWidths==L,1,'first'); % width when the end fo te freq. series is reached
77
78 %% 1st loop on filter width
79 for iiWidth = 1:lastWidth
80 NAverages = ceil( (iMaxForWidths(iiWidth)-iMin+1) / widths(iiWidth) ); % number of times the width avgWidth is going to be used
81 for kk=1:NAverages
82 %% second loop on filter iteration
83 iMax = min(L, iMin + widths(iiWidth)-1); % last sample for current average
84 idxAvg(iMin:iMax) = ones(1,iMax-iMin+1) * iAvg; % index of corresponding average
85 iMin = iMax + 1;
86 iAvg = iAvg + 1;
87 end
88 end
89
90 %% creating output data
91 sumMat = sparse(idxAvg, 1:L, ones(1,L));
92 nFreqs = sumMat*ones(L,1);
93
94 end