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)

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 linear 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.

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

Generated on Mon 03-Sep-2007 12:12:34 by m2html © 2003