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.m,v 1.4 2008/02/20 07:39:01 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.m,v 1.4 2008/02/20 07:39:01 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   else
0070     [y,coeffs] = ltpda_polyreg(x(xStart(i):xEnd(i)), order);
0071     dseg       = y.*win;
0072   end
0073   % compute fft
0074   xx   = fft(dseg,nfft);
0075   Sxx = ((Sxx .* (i-1)) +  xx.*conj(xx)) ./ i;
0076 end
0077 % take only one-sided
0078 Sxx = Sxx(1:nfft/2+1);
0079 
0080 % scale the data as required
0081 [Sxx, info] = scale_xx(Sxx, win, fs, scale, inunits);
0082 
0083 % grow frequency vector
0084 f = linspace(0, fs/2, length(Sxx));
0085 
0086 % complete info
0087 info.nsecs = Nx/fs;
0088 info.navs  = k;
0089 info.nolap = noverlap/fs;
0090 
0091 if nargout == 3
0092 
0093   varargout{1} = f;
0094   varargout{2} = Sxx;
0095   varargout{3} = info;
0096 
0097 
0098 else
0099   disp('usage: [f, pxx, info] = ltpda_psd(x, win, nolap, nfft, fs, scale, inunits, order)');
0100   error('# incorrect usage');
0101 end
0102 
0103 
0104 %-----------------------------------------------------------------------------------------------
0105 function [yy, info] = scale_xx(xx, win, fs, norm, inunits)
0106 
0107 nfft = length(win);
0108 S1   = sum(win);
0109 S2   = sum(win.^2);
0110 enbw = fs * S2 / (S1*S1);
0111 
0112 switch norm
0113   case 'ASD'
0114     yy = sqrt(xx) * sqrt(2 / (fs*S2));
0115     info.units = [sprintf('%s', inunits) '/\surdHz'];
0116   case 'PSD'
0117     yy = xx * 2 /(fs*S2);
0118     info.units = sprintf('%s^2/Hz', inunits);
0119   case 'AS'
0120     yy = sqrt(xx) * sqrt(2) / S1;
0121     info.units = sprintf('%s', inunits);
0122   case 'PS'
0123     yy = xx * 2 / (S1^2);
0124     info.units = sprintf('%s^2', inunits);
0125   case 'none'
0126     yy = xx;
0127     info.units = '';
0128   otherwise
0129     error('Unknown normalisation');
0130 end
0131 
0132 info.nfft = nfft;
0133 info.enbw = enbw;
0134 info.norm = norm;
0135 
0136 %-----------------------------------------------------------------------------------------------
0137 function [L,noverlap,win,msg] = segment_info(M,win,noverlap)
0138 %SEGMENT_INFO   Determine the information necessary to segment the input data.
0139 %
0140 %   Inputs:
0141 %      M        - An integer containing the length of the data to be segmented
0142 %      WIN      - A scalar or vector containing the length of the window or the window respectively
0143 %                 (Note that the length of the window determines the length of the segments)
0144 %      NOVERLAP - An integer containing the number of samples to overlap (may be empty)
0145 %
0146 %   Outputs:
0147 %      L        - An integer containing the length of the segments
0148 %      NOVERLAP - An integer containing the number of samples to overlap
0149 %      WIN      - A vector containing the window to be applied to each section
0150 %      MSG      - A string containing possible error messages
0151 %
0152 %
0153 %   The key to this function is the following equation:
0154 %
0155 %      K = (M-NOVERLAP)/(L-NOVERLAP)
0156 %
0157 %   where
0158 %
0159 %      K        - Number of segments
0160 %      M        - Length of the input data X
0161 %      NOVERLAP - Desired overlap
0162 %      L        - Length of the segments
0163 %
0164 %   The segmentation of X is based on the fact that we always know M and two of the set
0165 %   {K,NOVERLAP,L}, hence determining the unknown quantity is trivial from the above
0166 %   formula.
0167 
0168 % Initialize outputs
0169 L = [];
0170 msg = '';
0171 
0172 % Check that noverlap is a scalar
0173 if any(size(noverlap) > 1),
0174     msg.identifier = generatemsgid('invalidNoverlap');
0175     msg.message = 'You must specify an integer number of samples to overlap.';
0176     return;
0177 end
0178 
0179 if isempty(win),
0180     % Use 8 sections, determine their length
0181     if isempty(noverlap),
0182         % Use 50% overlap
0183         L = fix(M./4.5);
0184         noverlap = fix(0.5.*L);
0185     else
0186         L = fix((M+7.*noverlap)./8);
0187     end
0188     % Use a default window
0189     win = hamming(L);
0190 else
0191     % Determine the window and its length (equal to the length of the segments)
0192     if ~any(size(win) <= 1) | ischar(win),
0193         msg.identifier = generatemsgid('invalidWindow');
0194         msg.message = 'The WINDOW argument must be a vector or a scalar.';
0195         return
0196     elseif length(win) > 1,
0197         % WIN is a vector
0198         L = length(win);
0199     elseif length(win) == 1,
0200         L = win;
0201         win = hamming(win);
0202     end
0203     if isempty(noverlap),
0204         % Use 50% overlap
0205         noverlap = fix(0.5.*L);
0206     end
0207 end
0208 
0209 % Do some argument validation
0210 if L > M,
0211     msg.identifier = generatemsgid('invalidSegmentLength');
0212     msg.message = 'The length of the segments cannot be greater than the length of the input signal.';
0213     return;
0214 end
0215 
0216 if noverlap >= L,
0217     msg.identifier = generatemsgid('invalidNoverlap');
0218     msg.message = 'The number of samples to overlap must be less than the length of the segments.';
0219     return;
0220 end
0221 
0222 %-----------------------------------------------------------------------------------------------
0223 
0224 % END

Generated on Fri 07-Mar-2008 15:46:43 by m2html © 2003