Home > classes > @ao > welch.m

welch

PURPOSE ^

WELCH Welch spectral estimation method.

SYNOPSIS ^

function varargout = welch(varargin)

DESCRIPTION ^

 WELCH Welch spectral estimation method.

 [pxx, f, info] = welch(x,pl)
 [pxx, f, info] = welch(x,y,pl)

 PARAMETERS: 'Win'   - a specwin window object [default: Kaiser -200dB psll]
             'Olap' - segment percent overlap [default: taken from window function]
             'Nfft'  - number of samples in each fft [default: length of input data]
             'Scale' - one of
                                'ASD' - amplitude spectral density
                                'PSD' - power spectral density [default]
                                'AS'  - amplitude spectrum
                                'PS'  - power spectrum
                       * applies only to spectrum 'Type' 'psd'
             'Order' - order of segment detrending
                        -1 - no detrending
                         0 - subtract mean [default]
                         1 - subtract linear fit
                         N - subtract fit of polynomial, order N
             'Type'  - type of spectrum:
                       'psd'      - compute PSD
                       'cpsd'     - compute cross-spectral density
                       'tfe'      - estimate transfer function between inputs
                       'mscohere' - estimate cross-coherence

 Copied directly from MATLAB and extended to do segment-wise detrending
 and to take a plist input

 M Hewitson 08-05-08

 $Id: welch.m,v 1.5 2008/09/03 17:05:51 ingo Exp $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % WELCH Welch spectral estimation method.
0002 %
0003 % [pxx, f, info] = welch(x,pl)
0004 % [pxx, f, info] = welch(x,y,pl)
0005 %
0006 % PARAMETERS: 'Win'   - a specwin window object [default: Kaiser -200dB psll]
0007 %             'Olap' - segment percent overlap [default: taken from window function]
0008 %             'Nfft'  - number of samples in each fft [default: length of input data]
0009 %             'Scale' - one of
0010 %                                'ASD' - amplitude spectral density
0011 %                                'PSD' - power spectral density [default]
0012 %                                'AS'  - amplitude spectrum
0013 %                                'PS'  - power spectrum
0014 %                       * applies only to spectrum 'Type' 'psd'
0015 %             'Order' - order of segment detrending
0016 %                        -1 - no detrending
0017 %                         0 - subtract mean [default]
0018 %                         1 - subtract linear fit
0019 %                         N - subtract fit of polynomial, order N
0020 %             'Type'  - type of spectrum:
0021 %                       'psd'      - compute PSD
0022 %                       'cpsd'     - compute cross-spectral density
0023 %                       'tfe'      - estimate transfer function between inputs
0024 %                       'mscohere' - estimate cross-coherence
0025 %
0026 % Copied directly from MATLAB and extended to do segment-wise detrending
0027 % and to take a plist input
0028 %
0029 % M Hewitson 08-05-08
0030 %
0031 % $Id: welch.m,v 1.5 2008/09/03 17:05:51 ingo Exp $
0032 %
0033 
0034 %   Author(s): P. Pacheco
0035 %   Copyright 1988-2006 The MathWorks, Inc.
0036 %   $Revision: 1.5 $  $Date: 2008/09/03 17:05:51 $
0037 
0038 %   References:
0039 %     [1] Petre Stoica and Randolph Moses, Introduction To Spectral
0040 %         Analysis, Prentice-Hall, 1997, pg. 15
0041 %     [2] Monson Hayes, Statistical Digital Signal Processing and
0042 %         Modeling, John Wiley & Sons, 1996.
0043 
0044 % Compute the periodogram power spectrum of each segment and average always
0045 % compute the whole power spectrum, we force Fs = 1 to get a PS not a PSD.
0046 
0047 function varargout = welch(varargin)
0048 
0049   if nargin == 3
0050     a  = varargin{1};
0051     esttype = varargin{2};
0052     pl = varargin{3};
0053     x  = a.data.y;
0054     inunits = a.data.yunits;
0055   else
0056     a  = varargin{1};
0057     b  = varargin{2};
0058     esttype = varargin{3};
0059     pl = varargin{4};
0060     if a.data.fs ~= b.data.fs
0061       error('### Two time-series have different sample rates.');
0062     end
0063     inunits = b.data.yunits / a.data.yunits;
0064     x  = {a.data.y, b.data.y};
0065   end
0066 
0067   % Parse inputs
0068   win     = find(pl, 'Win');
0069   nfft    = find(pl, 'Nfft');
0070   olap    = find(pl, 'Olap')/100;
0071   scale   = find(pl, 'scale');
0072   Xolap   = floor(olap*nfft);
0073   fs      = a.data.fs;
0074   order   = find(pl, 'order');
0075 
0076   [x,M,isreal_x,y,Ly,win,winName,winParam,noverlap,k,L,options] = ...
0077     ao.welchparse(x,esttype,win.win, Xolap, nfft, fs);
0078 
0079   % Initialize
0080   Sxx = zeros(options.nfft,1);
0081 
0082   % Frequency vector was specified, return and plot two-sided PSD
0083   freqVectorSpecified = false; nrow = 1;
0084   if length(options.nfft) > 1,
0085     freqVectorSpecified = true;
0086     [ncol,nrow] = size(options.nfft);
0087   end
0088 
0089   % Compute the periodogram power spectrum of each segment and average always
0090   % compute the whole power spectrum, we force Fs = 1 to get a PS not a PSD.
0091 
0092   % Initialize
0093   if freqVectorSpecified,
0094     Sxx = zeros(length(options.nfft),1);
0095   else
0096     Sxx = zeros(options.nfft,1);
0097   end
0098   range = options.range;
0099 
0100   LminusOverlap = L-noverlap;
0101   xStart = 1:LminusOverlap:k*LminusOverlap;
0102   xEnd   = xStart+L-1;
0103   switch esttype,
0104     case {'ms','psd'},
0105       for i = 1:k,
0106         if order < 0
0107           seg = x(xStart(i):xEnd(i));
0108         else
0109           [seg,coeffs] = ltpda_polyreg(x(xStart(i):xEnd(i)), order);
0110         end
0111         [Sxxk,w] = ao.computeperiodogram(seg,win,...
0112           options.nfft,esttype,options.Fs);
0113         Sxx  = Sxx + Sxxk;
0114       end
0115     case 'cpsd',
0116       for i = 1:k,
0117         if order < 0
0118           xseg = x(xStart(i):xEnd(i));
0119         else
0120           [xseg,coeffs] = ltpda_polyreg(x(xStart(i):xEnd(i)), order);
0121         end
0122         if order < 0
0123           yseg = y(xStart(i):xEnd(i));
0124         else
0125           [yseg,coeffs] = ltpda_polyreg(y(xStart(i):xEnd(i)), order);
0126         end
0127         [Sxxk,w] =  ao.computeperiodogram({xseg,...
0128           yseg},win,options.nfft,esttype,options.Fs);
0129         Sxx  = Sxx + Sxxk;
0130       end
0131     case 'tfe'
0132       Sxy = zeros(options.nfft,1); % Initialize
0133       for i = 1:k,
0134         if order < 0
0135           xseg = x(xStart(i):xEnd(i));
0136         else
0137           [xseg,coeffs] = ltpda_polyreg(x(xStart(i):xEnd(i)), order);
0138         end
0139         if order < 0
0140           yseg = y(xStart(i):xEnd(i));
0141         else
0142           [yseg,coeffs] = ltpda_polyreg(y(xStart(i):xEnd(i)), order);
0143         end
0144         [Sxxk,w] = ao.computeperiodogram(xseg,...
0145           win,options.nfft,esttype,options.Fs);
0146         [Syxk,w] = ao.computeperiodogram({yseg,...
0147           x(xStart(i):xEnd(i))},win,options.nfft,esttype,options.Fs);
0148         Sxx  = Sxx + Sxxk;
0149         Sxy  = Sxy + Syxk;
0150       end
0151     case 'mscohere'
0152       % Note: (Sxy1+Sxy2)/(Sxx1+Sxx2) != (Sxy1/Sxy2) + (Sxx1/Sxx2)
0153       % ie, we can't push the computation of Cxy into computeperiodogram.
0154       Sxy = zeros(options.nfft,1); % Initialize
0155       Syy = zeros(options.nfft,1); % Initialize
0156       for i = 1:k,
0157         if order < 0
0158           xseg = x(xStart(i):xEnd(i));
0159         else
0160           [xseg,coeffs] = ltpda_polyreg(x(xStart(i):xEnd(i)), order);
0161         end
0162         if order < 0
0163           yseg = y(xStart(i):xEnd(i));
0164         else
0165           [yseg,coeffs] = ltpda_polyreg(y(xStart(i):xEnd(i)), order);
0166         end
0167         [Sxxk,w] = ao.computeperiodogram(xseg,...
0168           win,options.nfft,esttype,options.Fs);
0169         [Syyk,w] =  ao.computeperiodogram(yseg,...
0170           win,options.nfft,esttype,options.Fs);
0171         [Sxyk,w] = ao.computeperiodogram({xseg,...
0172           yseg},win,options.nfft,esttype,options.Fs);
0173         Sxx  = Sxx + Sxxk;
0174         Syy  = Syy + Syyk;
0175         Sxy  = Sxy + Sxyk;
0176       end
0177   end
0178 
0179   Sxx = Sxx./k; % Average the sum of the periodograms
0180   if any(strcmpi(esttype,{'tfe','mscohere'})),
0181     Sxy = Sxy./k; % Average the sum of the periodograms
0182 
0183     if strcmpi(esttype,'mscohere'),
0184       Syy = Syy./k; % Average the sum of the periodograms
0185     end
0186   end
0187 
0188   % Generate the freq vector directly in Hz to avoid roundoff errors due to
0189   % conversions later.
0190   if ~freqVectorSpecified,
0191     w = psdfreqvec('npts',options.nfft, 'Fs',options.Fs);
0192   else
0193     if strcmpi(options.range,'onesided')
0194       warning(generatemsgid('InconsistentRangeOption'),...
0195         'Ignoring ''onesided'' option. When a frequency vector is specified, a ''twosided'' PSD is computed');
0196     end
0197     options.range = 'twosided';
0198   end
0199 
0200 
0201   % Compute the 1-sided or 2-sided PSD [Power/freq] or mean-square [Power].
0202   % Also, corresponding freq vector and freq units.
0203   [Pxx,w,units] = computepsd(Sxx,w,options.range,options.nfft,options.Fs,esttype);
0204 
0205   if any(strcmpi(esttype,{'tfe','mscohere'})),
0206     % Cross PSD.  The frequency vector and xunits are not used.
0207     [Pxy,f,xunits] = computepsd(Sxy,w,options.range,options.nfft,options.Fs,esttype);
0208 
0209     % Transfer function estimate.
0210     if strcmpi(esttype,'tfe'),
0211       Pxx = Pxy ./ Pxx; % Txy
0212     end
0213 
0214     % Magnitude Square Coherence estimate.
0215     if strcmpi(esttype,'mscohere'),
0216       % Auto PSD for 2nd input vector. The freq vector & xunits are not
0217       % used.
0218       [Pyy,f,xunits] = computepsd(Syy,w,options.range,options.nfft,options.Fs,esttype);
0219       Pxx = (abs(Pxy).^2)./(Pxx.*Pyy); % Cxy
0220     end
0221   end
0222 
0223   % Scale to required units
0224   [Pxx, info] = ao.welchscale(Pxx, win, fs, scale, inunits);
0225   info.navs = k;
0226 
0227   varargout = {Pxx, w, info};
0228 end
0229

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003