前言

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

目錄

• 前言
• 目錄
• MTM實(shí)現(xiàn)
• 函數(shù)computepsd

MTM實(shí)現(xiàn)

檢查輸入 narginchk(1,13); len = length(varargin); 將輸入的x轉(zhuǎn)換為列向量 % convert vectors to column vectors for backwards compatibility if isvector(x)x = x(:); end 如果x的維度大于2，只支持向量和矩陣，則報(bào)錯(cuò) if numel(size(x))>2error(message('signal:pmtm:MustBeVectorOr2DMatrix')); end 是否丟掉最后一個(gè)窗 tf = strcmpi('droplasttaper',varargin); 指明第幾個(gè)輸入?yún)?shù)是丟窗 indx = find(tf==1); drop 的值只能是真或假 if (~isempty(indx) && ~islogical(varargin{indx+1}))error(message('signal:pmtm:MustBeLogical', 'Droplasttaper')); end 將這兩個(gè)參數(shù)，移到最后 % 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 解析輸入，這個(gè)函數(shù)在2種分析了 % Parse the inputs, set up default values, and return any error messages. params = parseinputs(x,varargin{:}); 強(qiáng)制使用精確規(guī)則 % 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計(jì)算功率譜 % Compute the two-sided power spectrum via MTM. [S,k,w] = mtm_spectrum(x,params); 生成頻率點(diǎn) % Generate the freq vector in correct units to avoid roundoff errors due to % converting frequency units later. nfft = params.nfft; [~,ncol] = size(nfft); 計(jì)算單邊或者雙邊功率譜或者均方譜 % Compute the 1-sided or 2-sided PSD [Power/freq] or mean-square [Power]. % Also, compute the corresponding freq vector & freq units. 第二個(gè)行之后的能量都*2/總頻率 赫茲數(shù) [Pxx,f,units] = computepsd(S,w,params.range,params.nfft,params.Fs,'psd'); 計(jì)算置信域，如果有需求的話。 % Calculate confidence limits ONLY when needed, since it can take a while. 強(qiáng)制雙進(jìn)度，在計(jì)算。 % Enforce double precision arithmetic on the calculations. chi2conf already 強(qiáng)制雙進(jìn)度計(jì)算 % 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 多少個(gè)信號(hào)源 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 如果沒有輸出參數(shù) 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結(jié)構(gòu)化的數(shù)據(jù)% Create a spectrum object to store in the PSD object's metadata. 存儲(chǔ)dsp數(shù)據(jù)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);頻率點(diǎn) 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和頻率點(diǎn)和置信度else% use legacy output ordervarargout = {Pxx,Pxxc,f};end end

函數(shù)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'??????????````

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。