前言

之前講了MTM（多錐形窗譜估計）的相關原理，現在來分析一下它的matlab實現。
想要復習的可以參考一下之前的文件：
現代譜估計：多窗口譜
想要復習一下如何實現的可以參考：
MTM：matlab實現1MTM：matlab實現1
MTM：matlab實現2參數解析MTM參數解析
MTM：matlab實現3譜功率計算MTM譜功率計算
MTM：matlab實現4主函數解析MTM 主函數

目錄

• 前言
• 目錄
• MTM實現
• 函數computepsd

MTM實現

檢查輸入 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))>2error(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 解析輸入，這個函數在2種分析了 % 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); 計算單邊或者雙邊功率譜或者均方譜 % Compute the 1-sided or 2-sided PSD [Power/freq] or mean-square [Power]. % Also, compute the corresponding freq vector & freq units. 第二個行之后的能量都*2/總頻率 赫茲數 [Pxx,f,units] = computepsd(S,w,params.range,params.nfft,params.Fs,'psd'); 計算置信域，如果有需求的話。 % Calculate confidence limits ONLY when needed, since it can take a while. 強制雙進度，在計算。 % Enforce double precision arithmetic on the calculations. chi2conf already 強制雙進度計算 % enforces double precision arithmetic. if ~strcmp(params.conflevel,'omitted') && (nargout==0 || nargout>2)% 'ConfidenceLevel' pv-pair specifiedNchan = 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 syntaxif ~strcmp(params.ConfInt,'omitted')c = chi2conf(params.ConfInt,k);else % use default value of .95k自由度卡方分布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); elsePxxc = []; 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 單位是hzhpsd = dspdata.psd([Pxx Pxxc],f{:},'SpectrumType',params.range); hpsd結構化的數據% Create a spectrum object to store in the PSD object's metadata. 存儲dsp數據hspec = spectrum.mtm(params.E,params.V,params.MTMethod);hpsd.Metadata.setsourcespectrum(hspec);if params.centerdccenterdc(hpsd);endhp = plot(hpsd);畫圖譜圖if 3*size(x,2)==numel(hp)nChan = size(x,2);color = get(hp,'Color');for i=1:nChanset(hp(nChan+2*i-1), 'Color',color{i}, 'LineStyle','-.');set(hp(nChan+2*i), 'Color',color{i}, 'LineStyle','-.');endend return end% center dc if needed if params.centerdc[Pxx, f, Pxxc] = psdcenterdc(Pxx, f, Pxxc, params); end中心化頻率 返回行向量作為輸出 % If the input is a vector and a row frequency vector was specified, % return output as a row vector for backwards compatibility. if ncol > 1 && nargout > 0 && isvector(x)Pxx = Pxx.';f = f.';% preserve (undocumented) behavior with legacy syntax.if strcmp(params.conflevel,'omitted') && nargout >= 3Pxxc = Pxxc.';end endif isa(Pxx,'single')% Cast to enforce precision rules.f = single(f);頻率點 end if nargout==1varargout = {Pxx};輸出pxx elseif nargout==2varargout = {Pxx,f};輸出pxx和f elseif nargout==3if ~strcmp(params.conflevel,'omitted')% use preferred output ordervarargout = {Pxx,f,Pxxc};輸出pxx和頻率點和置信度else% use legacy output ordervarargout = {Pxx,Pxxc,f};end end

函數computepsd

“`
Compute the one-sided or two-sided PSD or Mean-Square.
[Pxx,W,UNITS] = computepsd(Sxx,W,RANGE,NFFT,Fs,ESTTYPE) where the
inputs and outputs are:

Inputs:Sxx - Whole power spectrum [Power]; it can be a vector or a matrix.For matrices the operation is applied to each column.W - Frequency vector in rad/sample or in Hz.RANGE - Determines if a 'onesided' or a 'twosided' Pxx and Sxx arereturned.NFFT - Number of frequency points.Fs - Sampling Frequency.ESTTYPE - A string indicating the estimate type: 'psd', or 'ms' value.Outputs:Pxx - One-sided or two-sided PSD or MEAN-SQUARE (not scaled by Fs)depending on the input arguments RANGE and TYPE.W - Frequency vector 0 to 2*Nyquist or 0 to Nyquist depending onrange, units will be either rad/sample (if Fs is empty) or Hz(otherwise).UNITS -sample' or 'Hz'??????????````

總結

以上是生活随笔為你收集整理的MTM：matlab实现5主函数解码的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。