前言

之前講了MTM（多錐形窗譜估計）的相關原理，現在來分析一下它的matlab實現。
想要復習的可以參考一下之前的文件：
現代譜估計：多窗口譜

正文

首先是基本的子函數

findNFFT

從傳遞過來的參數里面，找到特定的NFFT或者頻率向量，返回給調用函數

------------------------------------------------------------------ function nfftTemp = findNFFT(varargin) % FINDNFFT Finds the specified NFFT or frequency vector from the optional % arguments passednfftTemp = []; for cnt = 1:length(varargin)if isnumeric(varargin{cnt}), nfftTemp = varargin{cnt};找到對應的參數后，立馬退出。break;end endend end

getEV

函數[E,V,NW,indx,nfft_temp,varargin]=getEV(N,varargin)
將輸入參數解析為 (但不包括) Nfft. 如果e和v未指定, 則計算 E 和 V

%---------------------------------------------------------------------- function [E,V,NW,indx,nfft_temp,varargin] = getEV(N,varargin) % GETEV Parse the input arguments up to, but not including, Nfft and % calculate E and V if not specified. % 輸入： N 輸入數據x的長度 % Inputs: % N - Length of the input data sequence, x.varargin -傳遞給pmtm的輸入參數列表，除了x% varargin - Input parameter list passed to pmtm, except for x. %輸出 % Outputs: e 離散扁球序列矩陣 % E - Matrix containing the discrete prolate spheroidal % sequences (dpss). v 包含dpss中心的向量% V - Vector containing the concentration of the dpss.NW 時間窗口是4，默認是4 % NW - Time-bandwidth product; default is 4.指示輸入參數列表開始選項的索引 % indx - Index indicating starting location of options in pmtm's % input argument list.特定的NFFT點數 % nfft_temp - NFFT or Frequency vector specified. Empty if not specified 設置默認的或者初始化輸出變量，避免早返回% Define defaults & initialize output variables (in case of early return). NW = 4; indx = 1; % Index where the options begin in the input arg list nfft_temp = []; 判斷是否輸入droplasttaper參數 tf = strcmpi('droplasttaper',varargin); 查找tf==1的位置。 loc = find(tf==1); 如果找到了droplasttaper， 默認是dlt=1，丟掉最后一幀。 if ~isempty(loc)dlt = varargin{loc+1}; % droplasttapervarargin = varargin(1:loc-1); elsedlt = true; % default value end第二種輸入參數給pmtm可以是以下的類型： 1（x,Nw,...）標量 x，矩陣e，向量v x，元胞數組包含dpss的所有參數 一個特定的功率譜選項% The 2nd input arg to pmtm can be a % 1. (X,NW,...) scalar % 2. (X,E,V,...) matrix E, hence, 3rd input must be a vector (V) % 3. (X,{dpss_params},...) cell containing the input argument list to dpss % 4. a string specifying a psdoption滿足下列條件，則NW=varargin『1』. if ~isempty(varargin) && ~ischar(varargin{1})if ~isempty(varargin{1})NW = varargin{1};endindx = 2;if iscell(NW), % NW is a cell array => dpss_paramsif (NW{1}<1.25 && dlt)error(message('signal:pmtm:insufficientTimebandwidthproduct', 'NW', '1.25', 'Droplasttaper', 'true'));end if (NW{1}<0.75 && ~dlt)error(message('signal:pmtm:insufficientTimebandwidthproduct', 'NW', '0.75', 'Droplasttaper', 'false')); end 生成對應的E,V值[E,V] = dpss(N,NW{:}); numvec = length(V);if dltif numvec > 2E = E(:,1:numvec-1);V = V(1:numvec-1);elseerror(message('signal:pmtm:inadequateNumtapers', '3', 'Droplasttaper', 'true'));endelseif numvec < 2error(message('signal:pmtm:inadequateNumtapers', '2', 'Droplasttaper', 'false'));endendNW值NW = NW{1};if nargin > 2, nfft_temp = findNFFT(varargin{2:end}); endE，v提前設計好。elseif length(NW)>1, % NW is the matrix E (==>V must be specified)E = NW;if length(varargin)<2,error(message('signal:pmtm:MustProvideVWithE', 'V', 'E'));elseV = varargin{2};if nargin > 3, nfft_temp = findNFFT(varargin{3:end}); endendnumvec = length(V);if size(E,2)~=numvecerror(message('signal:pmtm:MismatchEV', 'E', 'V'));end NW = size(E,2)/2; indx = 3; % Update index of beginning of options in the input arg list if dltif (numvec < 3)error(message('signal:pmtm:inadequateNumtapers', '3', 'Droplasttaper', 'true'));elseE = E(:,1:numvec-1);V = V(1:numvec-1);endelseif(numvec < 2)error(message('signal:pmtm:inadequateNumtapers', '2', 'Droplasttaper', 'false'));endend 如果NW是個標量else % NW is a scalarif (NW<1.25 && dlt)error(message('signal:pmtm:insufficientTimebandwidthproduct', 'NW', '1.25', 'Droplasttaper', 'true'));endif (NW<0.75 && ~dlt)error(message('signal:pmtm:insufficientTimebandwidthproduct', 'NW', '0.75', 'Droplasttaper', 'false'));end% Get the dpss, one way or another:[E,V] = dpss(N,NW);numvec = length(V);if dltif numvec > 2E = E(:,1:numvec-1);V = V(1:numvec-1);elseerror(message('signal:pmtm:inadequateNumtapers', '3', 'Droplasttaper', 'true'));endelseif numvec < 2error(message('signal:pmtm:inadequateNumtapers', '2', 'Droplasttaper', 'false'));endendif nargin > 2, nfft_temp = findNFFT(varargin{2:end}); endend else% Get the dpss, one way or another:[E,V] = dpss(N,NW);numvec = length(V);if dltif numvec > 2E = E(:,1:numvec-1);V = V(1:numvec-1);elseerror(message('signal:pmtm:inadequateNumtapers', '3', 'Droplasttaper', 'true'));endelseif numvec < 2error(message('signal:pmtm:inadequateNumtapers', '2', 'Droplasttaper', 'false'));endendnfft_temp = []; end

總結

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

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