多窗谱分析

leadersheng

narginchk(1,13);
len = length(varargin); %将输入的x转换为列向量
% convert vectors to column vectors for backwards compatibility
if isvector(x)
    x = x(:);
end
%如果x的维度大于2，只支持向量和矩阵，则报错
if numel(size(x))>2
    error(message('signal:pmtm:MustBeVectorOr2DMatrix'));
end
tf = strcmpi('droplasttaper',varargin);
%指明第几个输入参数是丢窗
indx = find(tf==1);
%drop 的值只能是真或假
if (~isempty(indx) && ~islogical(varargin{indx+1}))
    error(message('signal:pmtm:MustBeLogical', 'Droplasttaper'));
end
% If the 'droplasttaper' pv-pair is used, move it to the end of varargin
if (~isempty(indx) && (indx+1 ~= len))
    dummy = varargin(1:indx-1);
    dummy(indx:len-2) = varargin(indx+2:len);
    dummy(len-1:len) = varargin(indx:indx+1);
    varargin = dummy;
end
% Parse the inputs, set up default values, and return any error messages.
params = parseinputs(x,varargin{:});
%强制使用精确规则
% Cast to enforce Precision Rules
% Already checked for invalid character inputs (NW, NFFT,Fs) in 'parseinputs->psdoptions'
params.nfft = double(params.nfft);
params.Fs = double(params.Fs);
%使用mtm计算功率谱
% Compute the two-sided power spectrum via MTM.
[S,k,w] = mtm_spectrum(x,params);  
% Generate the freq vector in correct units to avoid roundoff errors due to
% converting frequency units later.
nfft = params.nfft;
[~,ncol] = size(nfft);
[Pxx,f,units] = computepsd(S,w,params.range,params.nfft,params.Fs,'psd');  
if ~strcmp(params.conflevel,'omitted') && (nargout==0 || nargout>2)
    Nchan = size(x,2);
    c = chi2conf(params.conflevel,k);
    Pxxc(:,2:2:2*Nchan) = double(Pxx)*c(2);
    Pxxc(:,1:2:2*Nchan) = double(Pxx)*c(1);
elseif (~strcmp(params.ConfInt,'omitted') && nargout==0) || nargout>2
    % using legacy syntax
    if ~strcmp(params.ConfInt,'omitted')
        c = chi2conf(params.ConfInt,k);
    else % use default value of 0.95
        c = chi2conf(.95,k);
    end
  Nchan = size(x,2);
  %置信度上下界
  Pxxc(:,2:2:2*Nchan) = double(Pxx)*c(2);
  Pxxc(:,1:2:2*Nchan) = double(Pxx)*c(1);
else
  Pxxc = [];
end
% Output 如果没有输出参数
if nargout==0
    %画出功率谱图
    % If no output arguments are specified plot the PSD w/ conf intervals.
    f = {f};
    if strcmpi(units,'Hz'), f = [f {'Fs',params.Fs}];
    end
    hpsd = dspdata.psd([Pxx Pxxc],f{:},'SpectrumType',params.range);
    hspec = spectrum.mtm(params.E,params.V,params.MTMethod);
    hpsd.Metadata.setsourcespectrum(hspec);

   if params.centerdc
     centerdc(hpsd);
   end

   hp = plot(hpsd);%画图谱图
   if 3*size(x,2)==numel(hp)
     nChan = size(x,2);
     color = get(hp,'Color');
     for i=1:nChan
       set(hp(nChan+2*i-1), 'Color',color{i}, 'LineStyle','-.');
       set(hp(nChan+2*i), 'Color',color{i}, 'LineStyle','-.');
     end
   end     
   return
end

% center dc if needed
if params.centerdc
   [Pxx, f, Pxxc] = psdcenterdc(Pxx, f, Pxxc, params);
end%%%%中心化频率
if ncol > 1 && nargout > 0 && isvector(x)
    Pxx = Pxx.';
    f = f.';
    % preserve (undocumented) behavior with legacy syntax.
    if strcmp(params.conflevel,'omitted') && nargout >= 3
        Pxxc = Pxxc.';
    end
end
if isa(Pxx,'single')
    % Cast to enforce precision rules.
    f = single(f);%频率点
end
if nargout==1
    varargout = {Pxx};
elseif nargout==2
    varargout = {Pxx,f};
elseif nargout==3
    if ~strcmp(params.conflevel,'omitted')
        % use preferred output order
        varargout = {Pxx,f,Pxxc};
    else % use legacy output order
        varargout = {Pxx,Pxxc,f};
    end
end