www.gusucode.com > signal 工具箱matlab源码程序 > signal/private/welch.m

    function varargout = welch(x,esttype,varargin)
%WELCH Welch spectral estimation method.
%   [Pxx,F] = WELCH(X,WINDOW,NOVERLAP,NFFT,Fs,SPECTRUMTYPE,TRACE,ESTTYPE)
%   [Pxx,F] = WELCH({X},WINDOW,NOVERLAP,NFFT,Fs,SPECTRUMTYPE,TRACE,'psd')
%   [Pxx,F] = WELCH({X},WINDOW,NOVERLAP,NFFT,Fs,SPECTRUMTYPE,TRACE,'ms')
%   [Pxy,F] = WELCH({X,Y},WINDOW,NOVERLAP,NFFT,Fs,SPECTRUMTYPE,'cpsd')
%   [Txy,F] = WELCH({X,Y},WINDOW,NOVERLAP,NFFT,Fs,SPECTRUMTYPE,'tfe')
%   [Cxy,F] = WELCH({X,Y},WINDOW,NOVERLAP,NFFT,Fs,SPECTRUMTYPE,'mscohere')
%   [Pxx,F,Pxxc] = WELCH(...)
%   [Pxx,F,Pxxc] = WELCH(...,'ConfidenceLevel',P)
%
%   Inputs:
%      see "help pwelch" for complete description of all input arguments.
%      ESTTYPE - is a string specifying the type of estimate to return, the
%                choices are: psd, cpsd, tfe, and mscohere.
%
%   Outputs:
%      Depends on the input string ESTTYPE:
%      Pxx - Power Spectral Density (PSD) estimate, or
%      MS  - Mean-square spectrum, or
%      Pxy - Cross Power Spectral Density (CPSD) estimate, or
%      Txy - Transfer Function Estimate (TFE), or
%      Cxy - Magnitude Squared Coherence.
%      F   - frequency vector, in Hz if Fs is specified, otherwise it has
%            units of rad/sample

%   Copyright 1988-2014 The MathWorks, Inc.

%   References:
%     [1] Petre Stoica and Randolph Moses, Introduction To Spectral
%         Analysis, Prentice-Hall, 1997, pg. 15
%     [2] Monson Hayes, Statistical Digital Signal Processing and
%         Modeling, John Wiley & Sons, 1996.

narginchk(2,10);
nargoutchk(0,3);

% Parse input arguments.
[x,~,~,y,~,win,winName,winParam,noverlap,k,L,options] = ...
    welchparse(x,esttype,varargin{:});
% Cast to enforce precision rules
options.nfft = signal.internal.sigcasttofloat(options.nfft,'double',...
  'WELCH','NFFT','allownumeric');
noverlap = signal.internal.sigcasttofloat(noverlap,'double','WELCH',...
  'NOVERLAP','allownumeric');
options.Fs = signal.internal.sigcasttofloat(options.Fs,'double','WELCH',...
  'Fs','allownumeric');
k = double(k);

if any([signal.internal.sigcheckfloattype(x,'single')...
    signal.internal.sigcheckfloattype(y,'single'),...
    isa(win,'single')])
  x = single(x);
  y = single(y);
  win = single(win);
end

% Frequency vector was specified, return and plot two-sided PSD
freqVectorSpecified = false; nrow = 1;
if length(options.nfft) > 1,
    freqVectorSpecified = true;
    [~,nrow] = size(options.nfft);
end

% Compute the periodogram power spectrum of each segment and average always
% compute the whole power spectrum, we force Fs = 1 to get a PS not a PSD.

% Initialize
if freqVectorSpecified,
    nFreqs = length(options.nfft);
    if strcmpi(options.range,'onesided')
        warning(message('signal:welch:InconsistentRangeOption'));
    end
    options.range = 'twosided';
else
    nFreqs = options.nfft;
end
% Cast to enforce precision rules
Sxx = zeros(nFreqs,size(x,2),class(x)); %#ok<*ZEROLIKE>

LminusOverlap = L-noverlap;
xStart = 1:LminusOverlap:k*LminusOverlap;
xEnd   = xStart+L-1;

switch esttype,
    case {'ms','power','psd'}
        % Combining method
        if options.maxhold,
            Sxx(:) = -Inf;
            cmethod = @(seg,nextseg) max(seg,k*nextseg); % k will be removed below
        elseif options.minhold,
            Sxx(:) = Inf;
            cmethod = @(seg,nextseg) min(seg,k*nextseg); % k will be removed below
        else
            cmethod = @plus;
        end
        [Pxx,w,units] = localComputeSpectra(Sxx,x,[],xStart,xEnd,win,...
            options,esttype,k,cmethod,freqVectorSpecified);
        
    case 'cpsd'
        numChan = max(size(x,2),size(y,2));
        Sxx = zeros(nFreqs,numChan,class(x)); % Initialize
        cmethod = @plus;
        [Pxx,w,units] = localComputeSpectra(Sxx,x,y,xStart,xEnd,win,...
            options,esttype,k,cmethod,freqVectorSpecified);
      
    case 'tfe'
        numChan = max(size(x,2),size(y,2));
        Sxy = zeros(nFreqs,numChan,class(x));
        cmethod = @plus;
        [Pxx,w,units] = localComputeSpectra(Sxx,x,[],xStart,xEnd,win,...
            options,esttype,k,cmethod,freqVectorSpecified);
        % Cross PSD.  The frequency vector and xunits are not used.
        Pxy = localComputeSpectra(Sxy,y,x,xStart,xEnd,win,...
            options,esttype,k,cmethod,freqVectorSpecified);

        Pxx = bsxfun(@rdivide, Pxy, Pxx);
        
    case 'mscohere'
        % Note: (Sxy1+Sxy2)/(Sxx1+Sxx2) != (Sxy1/Sxy2) + (Sxx1/Sxx2)
        % ie, we can't push the computation of Cxy into computeperiodogram.
        numChan = max(size(x,2),size(y,2));
        Sxy = zeros(nFreqs,numChan,class(x));
        Syy = zeros(nFreqs,size(y,2),class(x));
        cmethod = @plus;
        [Pxx,w,units] = localComputeSpectra(Sxx,x,[],xStart,xEnd,win,...
            options,esttype,k,cmethod,freqVectorSpecified);
        Pyy = localComputeSpectra(Syy,y,[],xStart,xEnd,win,...
            options,esttype,k,cmethod,freqVectorSpecified);
        % Cross PSD.  The frequency vector and xunits are not used.
        Pxy = localComputeSpectra(Sxy,x,y,xStart,xEnd,win,...
            options,esttype,k,cmethod,freqVectorSpecified);

        Pxx = (abs(Pxy).^2)./bsxfun(@times,Pxx,Pyy); % Cxy
end

% Compute confidence intervals if needed.
if ~strcmp(options.conflevel,'omitted')
    if any(strcmpi(esttype,{'ms','power','psd'})) && ~isequal(cmethod,@plus),
        % Always compute the confidence interval around the average spectrum
        cmethod = @plus;
        Sxxc = zeros(nFreqs,size(x,2),class(x));
        [avgPxx,w] = localComputeSpectra(Sxxc,x,[],xStart,xEnd,win,...
            options,esttype,k,cmethod,freqVectorSpecified);
    else
        avgPxx = Pxx;
    end
    Pxxc = confInterval(options.conflevel, avgPxx, x, w, options.Fs, k);
elseif nargout>2
    if any(strcmpi(esttype,{'ms','power','psd'})) && ~isequal(cmethod,@plus),
        % Always compute the confidence interval around the average spectrum
        cmethod = @plus;
        Sxxc = zeros(nFreqs,size(x,2),class(x));
        [avgPxx,w] = localComputeSpectra(Sxxc,x,[],xStart,xEnd,win,...
            options,esttype,k,cmethod,freqVectorSpecified);
    else
        avgPxx = Pxx;
    end
    Pxxc = confInterval(0.95, avgPxx, x, w, options.Fs, k);
else
    Pxxc = [];
end

if nargout==0
    w = {w};
    if strcmpi(units,'Hz'), w = [w,{'Fs',options.Fs}];  end
    % Create a spectrum object to store in the Data object's metadata.
    percOverlap = (noverlap/L)*100;
    hspec = spectrum.welch({winName,winParam},L,percOverlap); %#ok<DWELCH>
    
    switch lower(esttype)
        case 'tfe'
            if strcmpi(options.range,'onesided'), range='half'; else range='whole'; end
            h = dspdata.freqz(Pxx,w{:},'SpectrumRange',range);
        case 'mscohere'
            if strcmpi(options.range,'onesided'), range='half'; else range='whole'; end
            h = dspdata.magnitude(Pxx,w{:},'SpectrumRange',range);
        case 'cpsd'
            h = dspdata.cpsd(Pxx,w{:},'SpectrumType',options.range);
        case {'ms','power'}
            h = dspdata.msspectrum(Pxx,w{:},'SpectrumType',options.range); %#ok<DMSSPEC>
        otherwise
            h = dspdata.psd(Pxx,w{:},'SpectrumType',options.range); %#ok<DPPSD>
    end
    h.Metadata.setsourcespectrum(hspec);
    
    % plot the confidence levels if conflevel is specified.
    if ~isempty(Pxxc)
        h.ConfLevel = options.conflevel;
        h.ConfInterval = Pxxc;
    end
    % center dc component if specified
    if options.centerdc
        centerdc(h);
    end
    plot(h);
    if strcmp(esttype,'power')
        title(getString(message('signal:welch:WelchPowerSpectrumEstimate')));
    end
else
    if options.centerdc
        [Pxx, w, Pxxc] = psdcenterdc(Pxx, w, Pxxc, options);
    end
    
    % If the input is a vector and a row frequency vector was specified,
    % return output as a row vector for backwards compatibility.
    if nrow > 1 && isvector(x)
        Pxx = Pxx.'; w = w.';
    end
    
    % Cast to enforce precision rules   
    % Only cast if output is requested, otherwise, plot using double
    % precision frequency vector.
    if isa(Pxx,'single')
      w = single(w);
    end
    
    if isempty(Pxxc)
        varargout = {Pxx,w}; % Pxx=PSD, MEANSQUARE, CPSD, or TFE
    else
        varargout = {Pxx,w,Pxxc};
    end       
end

end

function [Pxx,w,units] = localComputeSpectra(Sxx,x,y,xStart,xEnd,win,options,esttype,k,cmethod,freqVectorSpecified)

if isempty(y),
    for ii = 1:k
        [Sxxk,w] = computeperiodogram(x(xStart(ii):xEnd(ii),:),win,...
            options.nfft,esttype,options.Fs);
        Sxx  = cmethod(Sxx,Sxxk);
    end
else
    for i = 1:k
        [Sxxk,w] =  computeperiodogram({x(xStart(i):xEnd(i),:),...
            y(xStart(i):xEnd(i),:)},win,options.nfft,esttype,options.Fs);
        Sxx  = cmethod(Sxx,Sxxk);
    end
end
Sxx = Sxx./k; % Average the sum of the periodograms

% Generate the freq vector directly in Hz to avoid roundoff errors due to
% conversions later.
if ~freqVectorSpecified
    w = psdfreqvec('npts',options.nfft, 'Fs',options.Fs);
end

% Compute the 1-sided or 2-sided PSD [Power/freq] or mean-square [Power].
% Also, corresponding freq vector and freq units.
[Pxx,w,units] = computepsd(Sxx,w,options.range,options.nfft,options.Fs,esttype);
end

function Pxxc = confInterval(CL, Pxx, x, w, fs, k)
%   Reference: D.G. Manolakis, V.K. Ingle and S.M. Kagon,
%   Statistical and Adaptive Signal Processing,
%   McGraw-Hill, 2000, Chapter 5
k = fix(k);
c = privatechi2conf(CL,k);

% Cast to enforce precision rules
Pxx = double(Pxx);

PxxcLower = Pxx*c(1);
PxxcUpper = Pxx*c(2);
Pxxc = reshape([PxxcLower; PxxcUpper],size(Pxx,1),2*size(Pxx,2));

% DC and Nyquist bins have only one degree of freedom for real signals
if isreal(x)
    realConf = chi2conf(CL,k/2);
    Pxxc(w == 0,1:2:end) = Pxx(w == 0,:) * realConf(1);
    Pxxc(w == 0,2:2:end) = Pxx(w == 0,:) * realConf(2);
    if isempty(fs)
        Pxxc(w == pi,1:2:end) = Pxx(w == pi,:) * realConf(1);
        Pxxc(w == pi,2:2:end) = Pxx(w == pi,:) * realConf(2);
    else
        Pxxc(w == fs/2,1:2:end) = Pxx(w == fs/2,:) * realConf(1);
        Pxxc(w == fs/2,2:2:end) = Pxx(w == fs/2,:) * realConf(2);
    end
end
end
% [EOF]