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.5 2008/03/31 07:03:12 mauro 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') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.5 2008/03/31 07:03:12 mauro 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