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') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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