前言

之前讲了MTM（多锥形窗谱估计）的相关原理，现在来分析一下它的matlab实现。
想要复习的可以参考一下之前的文件：
现代谱估计：多窗口谱
想要复习一下如何实现的可以参考：
MTM：matlab实现1MTM：matlab实现1
MTM：matlab实现2参数解析MTM参数解析

正文

子函数的讲解即将收尾，现在是最后一个，多窗口谱计算函数。
输入
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; it contains the 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 is 1.
%      range    - default is 'onesided' or real 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 of the dpss.
%      NW       - Time-bandwidth product; default is 4.
%
% 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 from the input structure for convenience.
nfft = params.nfft;
E  = params.E;
V  = params.V;
Fs = params.Fs;
如果没有采样频率，则默认为归一化频率2pi
if isempty(Fs)Fs = 2*pi;
end
列向量计算MTM，每一列是一个单独的采样。
行数即采样点数，列数即独立信号个数
N = size(x,1);
Nchan = size(x,2);
采用的多窗口的窗口个数
k = length(V);
如果nfft不是一个整数而是一个列向量，则代表它是想要计算的点。
if length(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 does not support matrices% use CZT to compute DFT on nfft evenly spaced samples around the% 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' after about 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 by the eigenvalues, as in Park et al. % Eqn. 9.; note that the eqn. apparently has a typo; as the weights% should be V and not 1/V.if strcmp(params.MTMethod,'eigen')wt = V(:);    % Park estimateelsewt = ones(k,1);endSchan = Sk*wt/k;endS(:,chan) = Schan;
end

MTM：matlab实现3谱功率计算相关推荐

ar谱matlab实验,用MATLAB进行AR模型功率谱分析
用MATLAB 进行AR 模型功率谱分析随机信号序列x(n)是均值为0方差为1的高斯型白噪声经过AR 模型 ()4 3219606.01697.29403.22137.211 ----+-+-= z ...
用matlab求解信号的DFT,利用MATLAB实现信号DFT的计算
07级电信(2)班刘坤洋 24 实验一利用MATLAB 实现信号DFT 的计算一.实验目的: 1.熟悉利用MATLAB 计算信号DFT 的方法 2.掌握利用MATLAB 实现由DFT 计算线性卷 ...
matlab 简单算例,(简单算例)基于Matlab的电力系统潮流编程计算.pdf
(简单算例)基于Matlab的电力系统潮流编程计算基于Matlab的电力系统潮流编程计算口黄扬威吴喜春郭志峰张斯翔 (三峡大学电气与新能源学院湖北·宜昌443002) 摘要:通过介绍电力系统的实际 ...
电气潮流运算Matlab怎么编程,基于Matlab的电力系统潮流编程计算
计算技术信息发展与 64 -- 科协论坛 · 2011 年第 6 期(下) -- 基于 Matlab 的电力系统潮流编程计算 □ 黄扬威吴喜春郭志峰张斯翔 (三峡大学电气与新能源学院湖北· ...
基于matlab的自动识别谱峰的程序设计,基于matlab的自动识别谱峰的程序设计毕业论文-资源下载人人文库网...
基于matlab的自动识别谱峰的程序设计毕业论文目录摘要1一绪论211几种常用寻峰方法的简单说明212小波变换413MATLAB小波分析工具箱6二小波分析基本原理721一维连续小波分析722一维离 ...
matlab ar谱分析,用MATLAB进行AR模型功率谱分析
用MATLAB进行AR模型功率谱分析用MATLAB进行AR模型功率谱分析随机信号序列x(n)是均值为0方差为1的高斯型白噪声经过AR模型 H z 1 1 2.2137z 1 2.9403z 2 2 ...
根据功率计算用电量和电费
欢迎大家关注笔者,你的关注是我持续更博的最大动力原创文章,转载告知,盗版必究根据功率计算用电量和电费文章目录: 1 度的概念 2 根据设备功率计算耗电度数 3 计算电费 1 度的概念功率=电流 ...
基于matlab的图解粒度参数计算,基于MATLAB的图解粒度参数计算-热带地理.PDF
基于MATLAB的图解粒度参数计算-热带地理第 26卷第 3期热带地理 Vol26,No3 2006年 8月 TROP ICAL GEO GRA PHY Aug. , 2006 基于 MA ...
如何根据灰度直方图计算标准差_如何根据电器功率计算电线的粗细？
一般来说,测算电线的粗细,需要根据功率计算电流,根据电流选择导线截面,根据导线的截面,导线或电缆的型号查厂家的该型号的导线电缆的直径.这里就涉及了:电线粗细与功率之间的关系计算:导线截面积与载流量的计 ...

MTM：matlab实现3谱功率计算

前言

目录

正文

MTM：matlab实现3谱功率计算相关推荐

最新文章

热门文章