Home > m > sigproc > frequency_domain > ltpda_psd.m

ltpda_psd

PURPOSE ^

LTPDA_PSD a psd estimation based on pwelch from MATLAB.

SYNOPSIS ^

function varargout = ltpda_psd(x, win, nolap, nfft, fs, scale, inunits, order)

DESCRIPTION ^

 LTPDA_PSD a psd estimation based on pwelch from MATLAB.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: LTPDA_PSD a psd estimation based on pwelch from MATLAB.
              The difference is that detrending is done on each
              windowed data segment.

 CALL:       [f, pxx, info] = ltpda_psd(x, win, nolap, nfft, fs, scale, inunits)

 INPUTS:     x      - data vector
             win    - a window function
                      Specify a vector of window samples or a length. If a length
                      is specified, a Hamming window is used.
             nfft   - number of samples per fft
             nolap  - number of samples to overlap each segment by
             fs     - sample rate of input time-series
             scale  - specify one of 'ASD', 'AS', 'PSD', 'PS'
           inunits  - Amplitude units of input time-series data.
            order   - order of detrending

 OUTPUTS:    f      - a frequency vector
             Sxx    - spectral estimate in units specified by 'scale'
             info   - a structure of information

 NOTE:       segment_info() is copied directly from MATLAB's welshparse.m.

 VERSION:    $Id: ltpda_psd.html,v 1.14 2008/03/31 10:27:39 hewitson Exp $

 HISTORY:    25-09-2006 M Hewitson
                Creation

 The following call returns a parameter list object that contains the
 default parameter values:

 >> pl = ltpda_psd('Params')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_psd(x, win, nolap, nfft, fs, scale, inunits, order)
0002 % LTPDA_PSD a psd estimation based on pwelch from MATLAB.
0003 %
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %
0006 % DESCRIPTION: LTPDA_PSD a psd estimation based on pwelch from MATLAB.
0007 %              The difference is that detrending is done on each
0008 %              windowed data segment.
0009 %
0010 % CALL:       [f, pxx, info] = ltpda_psd(x, win, nolap, nfft, fs, scale, inunits)
0011 %
0012 % INPUTS:     x      - data vector
0013 %             win    - a window function
0014 %                      Specify a vector of window samples or a length. If a length
0015 %                      is specified, a Hamming window is used.
0016 %             nfft   - number of samples per fft
0017 %             nolap  - number of samples to overlap each segment by
0018 %             fs     - sample rate of input time-series
0019 %             scale  - specify one of 'ASD', 'AS', 'PSD', 'PS'
0020 %           inunits  - Amplitude units of input time-series data.
0021 %            order   - order of detrending
0022 %
0023 % OUTPUTS:    f      - a frequency vector
0024 %             Sxx    - spectral estimate in units specified by 'scale'
0025 %             info   - a structure of information
0026 %
0027 % NOTE:       segment_info() is copied directly from MATLAB's welshparse.m.
0028 %
0029 % VERSION:    $Id: ltpda_psd.html,v 1.14 2008/03/31 10:27:39 hewitson Exp $
0030 %
0031 % HISTORY:    25-09-2006 M Hewitson
0032 %                Creation
0033 %
0034 % The following call returns a parameter list object that contains the
0035 % default parameter values:
0036 %
0037 % >> pl = ltpda_psd('Params')
0038 %
0039 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0040 
0041 % get segment information
0042 Nx = length(x);
0043 [L,noverlap,win,msg] = segment_info(Nx,win,nolap);
0044 % Compute the number of segments
0045 k = (Nx-noverlap)./(L-noverlap);
0046 
0047 % Uncomment the following line to produce a warning each time the data
0048 % segmentation does not produce an integer number of segements.
0049 if fix(k) ~= k
0050    warning('The number of segments is not an integer, truncating data.');
0051 end
0052 k = fix(k);
0053 
0054 % fix vectors
0055 x   = x(:);
0056 win = win(:);
0057 
0058 % Initialise the spectral estimates
0059 Sxx = zeros(nfft,1);
0060 
0061 % loop over segments
0062 LminusOverlap = L-noverlap;
0063 xStart = 1:LminusOverlap:k*LminusOverlap;
0064 xEnd   = xStart+L-1;
0065 
0066 for i = 1:k,
0067   % data segment, detrend and window
0068   if order < 0
0069     y       = x(xStart(i):xEnd(i)).*win;
0070   else
0071     [y,coeffs] = ltpda_polyreg(x(xStart(i):xEnd(i)), order);    
0072   end
0073   dseg       = y.*win;
0074   % compute fft
0075   xx   = fft(dseg,nfft);
0076   Sxx = ((Sxx .* (i-1)) +  xx.*conj(xx)) ./ i;
0077 end
0078 % take only one-sided
0079 Sxx = Sxx(1:nfft/2+1);
0080 
0081 % scale the data as required
0082 [Sxx, info] = scale_xx(Sxx, win, fs, scale, inunits);
0083 
0084 % grow frequency vector
0085 f = linspace(0, fs/2, length(Sxx));
0086 
0087 % complete info
0088 info.nsecs = Nx/fs;
0089 info.navs  = k;
0090 info.nolap = noverlap/fs;
0091 
0092 if nargout == 3
0093 
0094   varargout{1} = f;
0095   varargout{2} = Sxx;
0096   varargout{3} = info;
0097 
0098 
0099 else
0100   disp('usage: [f, pxx, info] = ltpda_psd(x, win, nolap, nfft, fs, scale, inunits, order)');
0101   error('# incorrect usage');
0102 end
0103 
0104 
0105 %-----------------------------------------------------------------------------------------------
0106 function [yy, info] = scale_xx(xx, win, fs, norm, inunits)
0107 
0108 nfft = length(win);
0109 S1   = sum(win);
0110 S2   = sum(win.^2);
0111 enbw = fs * S2 / (S1*S1);
0112 
0113 switch norm
0114   case 'ASD'
0115     yy = sqrt(xx) * sqrt(2 / (fs*S2));
0116     info.units = [sprintf('%s', inunits) '/\surdHz'];
0117   case 'PSD'
0118     yy = xx * 2 /(fs*S2);
0119     info.units = sprintf('%s^2/Hz', inunits);
0120   case 'AS'
0121     yy = sqrt(xx) * sqrt(2) / S1;
0122     info.units = sprintf('%s', inunits);
0123   case 'PS'
0124     yy = xx * 2 / (S1^2);
0125     info.units = sprintf('%s^2', inunits);
0126   case 'none'
0127     yy = xx;
0128     info.units = '';
0129   otherwise
0130     error('Unknown normalisation');
0131 end
0132 
0133 info.nfft = nfft;
0134 info.enbw = enbw;
0135 info.norm = norm;
0136 
0137 %-----------------------------------------------------------------------------------------------
0138 function [L,noverlap,win,msg] = segment_info(M,win,noverlap)
0139 %SEGMENT_INFO   Determine the information necessary to segment the input data.
0140 %
0141 %   Inputs:
0142 %      M        - An integer containing the length of the data to be segmented
0143 %      WIN      - A scalar or vector containing the length of the window or the window respectively
0144 %                 (Note that the length of the window determines the length of the segments)
0145 %      NOVERLAP - An integer containing the number of samples to overlap (may be empty)
0146 %
0147 %   Outputs:
0148 %      L        - An integer containing the length of the segments
0149 %      NOVERLAP - An integer containing the number of samples to overlap
0150 %      WIN      - A vector containing the window to be applied to each section
0151 %      MSG      - A string containing possible error messages
0152 %
0153 %
0154 %   The key to this function is the following equation:
0155 %
0156 %      K = (M-NOVERLAP)/(L-NOVERLAP)
0157 %
0158 %   where
0159 %
0160 %      K        - Number of segments
0161 %      M        - Length of the input data X
0162 %      NOVERLAP - Desired overlap
0163 %      L        - Length of the segments
0164 %
0165 %   The segmentation of X is based on the fact that we always know M and two of the set
0166 %   {K,NOVERLAP,L}, hence determining the unknown quantity is trivial from the above
0167 %   formula.
0168 
0169 % Initialize outputs
0170 L = [];
0171 msg = '';
0172 
0173 % Check that noverlap is a scalar
0174 if any(size(noverlap) > 1),
0175     msg.identifier = generatemsgid('invalidNoverlap');
0176     msg.message = 'You must specify an integer number of samples to overlap.';
0177     return;
0178 end
0179 
0180 if isempty(win),
0181     % Use 8 sections, determine their length
0182     if isempty(noverlap),
0183         % Use 50% overlap
0184         L = fix(M./4.5);
0185         noverlap = fix(0.5.*L);
0186     else
0187         L = fix((M+7.*noverlap)./8);
0188     end
0189     % Use a default window
0190     win = hamming(L);
0191 else
0192     % Determine the window and its length (equal to the length of the segments)
0193     if ~any(size(win) <= 1) | ischar(win),
0194         msg.identifier = generatemsgid('invalidWindow');
0195         msg.message = 'The WINDOW argument must be a vector or a scalar.';
0196         return
0197     elseif length(win) > 1,
0198         % WIN is a vector
0199         L = length(win);
0200     elseif length(win) == 1,
0201         L = win;
0202         win = hamming(win);
0203     end
0204     if isempty(noverlap),
0205         % Use 50% overlap
0206         noverlap = fix(0.5.*L);
0207     end
0208 end
0209 
0210 % Do some argument validation
0211 if L > M,
0212     msg.identifier = generatemsgid('invalidSegmentLength');
0213     msg.message = 'The length of the segments cannot be greater than the length of the input signal.';
0214     return;
0215 end
0216 
0217 if noverlap >= L,
0218     msg.identifier = generatemsgid('invalidNoverlap');
0219     msg.message = 'The number of samples to overlap must be less than the length of the segments.';
0220     return;
0221 end
0222 
0223 %-----------------------------------------------------------------------------------------------
0224 
0225 % END

Generated on Mon 31-Mar-2008 12:20:24 by m2html © 2003