子函數的講解即將收尾，現在是最后一個，多窗口譜計算函數。 輸入 x 數據向量 params，pmtm傳遞的參數，除了輸入x 它包含以下值域 nfft，評估功率譜的具體頻率點 Fs 采樣頻率 range 單邊還是雙邊 ConfInt 置信間隔 默認0.95 mtmmethod多窗口譜方法，默認自適應 E，dpss矩陣 V 包含dpss序列中心的向量 NW 時間帶寬乘積 輸出 s MTM計算出的功率譜 k 構建pxx的加窗數 w dft計算的頻率點
%----------------------------------------------------------------------
function [S,k,w] = mtm_spectrum(x,params)
%MTM_SPECTRUM Compute the power spectrum via MTM.
%
% Inputs:
% x - Input data vector.
% params - Structure containing pmtm's input parameter list, except for
% the input data sequence, x; itcontainsthe following fields:
% nfft - Number of frequency points to evaluate the PSD at;
% the default is max(256,2^nextpow2(N)).
% Fs - The sampling frequency; default is1.
% range - default is 'onesided' orreal signals and 'twosided' for
% - complex signals.
% ConfInt - Confidence interval; default is.95.
% MTMethod - Algorithm used in MTM; default is 'adapt'.
% E - Matrix containing the discrete prolate spheroidal
% sequences (dpss).
% V - Vector containing the concentration ofthe dpss.
% NW - Time-bandwidth product; default is4.
%
% Outputs:
% S - Power spectrum computed via MTM.
% k - Number of sequences used to form Pxx.
% w - Frequency vector for which DFT is calculated從輸入中提取相關的參數
% Extract some parameters fromthe input structure for convenience.
nfft = params.nfft;
E = params.E;
V = params.V;
Fs = params.Fs;
如果沒有采樣頻率，則默認為歸一化頻率2piif isempty(Fs)Fs = 2*pi;
end
列向量計算MTM，每一列是一個單獨的采樣。
行數即采樣點數，列數即獨立信號個數
N = size(x,1);
Nchan = size(x,2);
采用的多窗口的窗口個數
k = length(V);
如果nfft不是一個整數而是一個列向量，則代表它是想要計算的點。
iflength(nfft) > 1, isfreqVector = true; nfft_mod = length(nfft);
else isfreqVector = false;nfft_mod = nfft;
end
如果x是single，則信號功率譜也是single
if isa(x,'single') || isa(E,'single')S = zeros(nfft_mod, Nchan, 'single');
elseS = zeros(nfft_mod, Nchan);
end
對于每一個信號源
for chan=1:Nchan
計算加窗 離散傅里葉% Compute the windowed DFTs.如果指定輸入頻率向量或者不指定且N采樣個數小于設定的NFFT點數if (~isfreqVector && N<=nfft) || isfreqVector 計算每一列的加窗fft變換% Compute DFT using FFT or Goertzel使用基本函數，快速計算加窗后的離散傅里葉變換[Xx,w] = computeDFT(bsxfun(@times,E(:,1:k),x(:,chan)),nfft(:),Fs); sk是xx的能量值 Sk = abs(Xx).^2;else % Wrap the data modulo nfft if N > nfft. Note we cannot use datawrap % and FFT because datawrap doesnot support matrices% use CZT to compute DFT on nfft evenly spaced samples aroundthe% unit circle:Sk = abs(czt(bsxfun(@times,E(:,1:k),x(:,chan)),nfft(:))).^2;
獲得對應的頻率點w = psdfreqvec('npts',nfft,'Fs',Fs);end
計算多窗口譜估計，在0：NFFt上計算整個特征譜。% Compute the MTM spectral estimates, compute the whole spectrum 0:nfft.switch params.MTMethod,
自適應的情況：case 'adapt'
設置決定自適應的權重的代數% Set up the iteration to determine the adaptive weights: 初始功率譜，第一列的平方和/數量sig2=x(:,chan)'*x(:,chan)/N; % Power加窗1和加窗2的平均值作為初始值Schan=(Sk(:,1)+Sk(:,2))/2; % Initial spectrum estimate S1=zeros(nfft_mod,1); % The algorithm converges so fast that results are% usually 'indistinguishable' afterabout three iterations.% This version uses the equations from [2] (P&W pp 368-370).% Set tolerance for acceptance of spectral estimate:tol=.0005*sig2/nfft_mod;a=bsxfun(@times,sig2,(1-V));% Do the iteration:while sum(abs(Schan-S1)/nfft_mod)>tol% calculate weightsb=(Schan*ones(1,k))./(Schan*V'+ones(nfft_mod,1)*a'); % calculate new spectral estimatewk=(b.^2).*(ones(nfft_mod,1)*V');S1=sum(wk'.*Sk')./ sum(wk,2)';S1=S1';Stemp=S1; S1=Schan; Schan=Stemp; % swap S and S1endcase {'unity','eigen'}% Compute the averaged estimate: simple arithmetic averaging is used. % The Sk can also be weighted bythe eigenvalues, asin Park et al. % Eqn. 9.; note thatthe eqn. apparently has a typo; asthe weights% should be V andnot1/V.if strcmp(params.MTMethod,'eigen')wt = V(:); % Park estimateelsewt = ones(k,1);endSchan = Sk*wt/k;endS(:,chan) = Schan;
end
