--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/m/sigproc/frequency_domain/ltpda_spsd.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,94 @@
+% ltpda_spsd smooths a spectrum.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: ltpda_spsd smooths a spectrum.
+%
+% CALL:        [freqsOut, pow, nFreqs, sigmaP, sumMat] = ltpda_spsd(freqs, pow, linCoef, logCoef, sumMat,nFreqs)
+%
+% INPUTS:      freqsOut - frequency vector, can be left empty
+%              pow     - power spectrum (density)
+%              nFreqs  - frequency intervals used to derive averages
+%              linCoef, logCoef
+%                      - values to use to scale the smoothing
+%                        averaging filter. It will be in linCoef.NBins^logCoef.
+%                        logCoef should be 2/3 for PSD and 4/5 for PSD data
+%                        which is later plotted in ASD scale.
+%              mode   - 'keepFrequencies' convoles using a filter, leaving all
+%                        frequencies but they are correlated
+%                       'removeFrequencies' convoles using a filter, leaving
+%                        only data at some uncorrelated frequencies
+%                       'addUp' like above, but sums up instead of takin a
+%                        mean value.
+%
+% OUTPUTS:     freqs, pow : output frequency and power vector.
+%
+% <a href="matlab:web(ao.getInfo('spsd').tohtml, '-helpbrowser')">Parameter Sets</a>
+%
+% VERSION:     $Id: ltpda_spsd.m,v 1.4 2011/02/21 16:55:21 adrien Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function varargout = ltpda_spsd(varargin)
+  %% collecting inputs
+  freqs = varargin{1};
+  pow = varargin{2};
+  L = length(pow);
+  pow = reshape(pow, [L, 1]);
+  freqs = reshape(freqs, [numel(freqs), 1]);
+  linCoef = varargin{3};
+  logCoef = varargin{4};
+  if nargin >4
+    % optional inputs include the averaging matrix
+    sumMat = varargin{5};
+    nFreqs = varargin{6};
+  else
+    % averaging matrix must be computed from input data
+    [nFreqs, sumMat] = getMatrix(L, linCoef, logCoef);
+  end
+  %% evaluating output
+  computeVar = nargout>=4;
+  computeFreqs = ~isempty(freqs);
+  sumPow = sumMat * pow; % power sums
+  powOut = sumPow ./ nFreqs; % power average
+  if computeFreqs
+    freqsOut = (sumMat * freqs) ./ nFreqs; % correponding mean frequency
+  else
+    freqsOut = [];
+  end
+  if computeVar
+    nBinsEff = sumPow.^2 ./ (sumMat * pow.^2) ; % number of Chi2 variables for STD
+  else
+    nBinsEff = [];
+  end
+  
+  %% allocating output
+  varargout = {freqsOut, powOut, nFreqs, nBinsEff, sumMat};
+  
+end
+
+
+function [nFreqs, sumMat] = getMatrix(L, linCoef, logCoef)
+  %% initializing loop
+  iAvg = 1;
+  iMin = 1;
+  idxAvg = zeros(1,L); % index of corresponding average
+  widths = ceil(linCoef):ceil(linCoef * L^logCoef + 2); % proposed width of averagin filter
+  iMaxForWidths = min(L, floor( (widths./linCoef).^(1/logCoef) + widths-1 ) ); % last sample to be processed with a given filter width
+  lastWidth = find(iMaxForWidths==L,1,'first'); % width when the end fo te freq. series is reached
+  
+  %% 1st loop on filter width
+  for iiWidth = 1:lastWidth
+    NAverages = ceil( (iMaxForWidths(iiWidth)-iMin+1) / widths(iiWidth) ); % number of times the width avgWidth is going to be used
+    for kk=1:NAverages
+      %% second loop on filter iteration
+      iMax = min(L, iMin + widths(iiWidth)-1); % last sample for current average
+      idxAvg(iMin:iMax) = ones(1,iMax-iMin+1) * iAvg; % index of corresponding average
+      iMin = iMax + 1;
+      iAvg = iAvg + 1;
+    end
+  end
+
+  %% creating output data  
+  sumMat = sparse(idxAvg, 1:L, ones(1,L));
+  nFreqs = sumMat*ones(L,1);
+  
+end
\ No newline at end of file