子函數(shù)的講解即將收尾，現(xiàn)在是最后一個(gè)，多窗口譜計(jì)算函數(shù)。 輸入 x 數(shù)據(jù)向量 params，pmtm傳遞的參數(shù)，除了輸入x 它包含以下值域 nfft，評估功率譜的具體頻率點(diǎn) Fs 采樣頻率 range 單邊還是雙邊 ConfInt 置信間隔 默認(rèn)0.95 mtmmethod多窗口譜方法，默認(rèn)自適應(yīng) E，dpss矩陣 V 包含dpss序列中心的向量 NW 時(shí)間帶寬乘積 輸出 s MTM計(jì)算出的功率譜 k 構(gòu)建pxx的加窗數(shù) w dft計(jì)算的頻率點(diǎn)
%----------------------------------------------------------------------
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從輸入中提取相關(guān)的參數(shù)
% Extract some parameters fromthe input structure for convenience.
nfft = params.nfft;
E = params.E;
V = params.V;
Fs = params.Fs;
如果沒有采樣頻率，則默認(rèn)為歸一化頻率2piif isempty(Fs)Fs = 2*pi;
end
列向量計(jì)算MTM，每一列是一個(gè)單獨(dú)的采樣。
行數(shù)即采樣點(diǎn)數(shù)，列數(shù)即獨(dú)立信號個(gè)數(shù)
N = size(x,1);
Nchan = size(x,2);
采用的多窗口的窗口個(gè)數(shù)
k = length(V);
如果nfft不是一個(gè)整數(shù)而是一個(gè)列向量，則代表它是想要計(jì)算的點(diǎn)。
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
對于每一個(gè)信號源
for chan=1:Nchan
計(jì)算加窗 離散傅里葉% Compute the windowed DFTs.如果指定輸入頻率向量或者不指定且N采樣個(gè)數(shù)小于設(shè)定的NFFT點(diǎn)數(shù)if (~isfreqVector && N<=nfft) || isfreqVector 計(jì)算每一列的加窗fft變換% Compute DFT using FFT or Goertzel使用基本函數(shù)，快速計(jì)算加窗后的離散傅里葉變換[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;
獲得對應(yīng)的頻率點(diǎn)w = psdfreqvec('npts',nfft,'Fs',Fs);end
計(jì)算多窗口譜估計(jì)，在0：NFFt上計(jì)算整個(gè)特征譜。% Compute the MTM spectral estimates, compute the whole spectrum 0:nfft.switch params.MTMethod,
自適應(yīng)的情況：case 'adapt'
設(shè)置決定自適應(yīng)的權(quán)重的代數(shù)% Set up the iteration to determine the adaptive weights: 初始功率譜，第一列的平方和/數(shù)量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
