如何矢量化多锥度谱的计算？答案

【问题标题】：How can I vectorize the calculation of a multi-taper spectrum?如何矢量化多锥度谱的计算？
【发布时间】：2013-03-03 21:44:11
【问题描述】：

通常我使用 pmtm 计算信号的频谱：

signal = rand(1000,1);
NW = 4;
Fr = 1:50;
Fs = 200;
[p, fr] = pmtm( signal, NW, Fr, Fs);

但是，我正在寻找一种方法来对其进行矢量化处理，以便可以同时计算多个光谱。我试过了：

signal = rand(1000,10); %<--- notice I have 10 columns instead of 1
NW = 4;
Fr = 1:50;
Fs = 200;
[p, fr] = pmtm( signal, NW, Fr, Fs);

但它产生的错误并不能真正告诉我我做错了什么。我知道我可以循环调用pmtm。

这是错误：

使用错误。* 矩阵尺寸必须一致。

pmtm 中的错误>mtm_spectrum（第 231 行） [Xx,w] = computeDFT(E(:,1:k).*x(:,ones(1,k)),nfft,Fs);

pmtm 中的错误（第 142 行）[S,k,w] = mtm_spectrum(x,params);

这让我怀疑没有一种矢量化的方式来实现我想要的。我希望这里有人知道如何做到这一点。

【问题讨论】：

mtm_spectrum 强制x 成为x = x(:); 行中的向量。如果x 是矩阵，它将有 n_col 更大的维度，这就是你得到矩阵维度错误的原因。你试过arrayfun吗？它是伪矢量化解决方案，但可能会稍微加快执行时间。像arrayfun(@(x)pmtm(signal(:,x), NW, Fr, Fs),1:size(signal,2)) 这样的东西。但是你只能得到一个输出变量。
看起来E 和x 的尺寸不完全相同...
你可能无法避免这里的 for 循环，因为pmtm 很慢，所以这非常昂贵。如果您关心性能，那么您可以通过使用 dpss 预先计算 slepian 序列并将它们传递给 pmtm 来显着加快速度。

标签： matlab fft spectrum

【解决方案1】：

是什么让您认为 pmtm 是为处理矩阵而设计的？帮助菜单特别需要一个向量：

>> help pmtm
 pmtm   Power Spectral Density (PSD) estimate via the Thomson multitaper 
    method (MTM).
    Pxx = pmtm(X) returns the PSD of a discrete-time signal vector X in 
    the vector Pxx.  Pxx is the distribution of power per unit frequency.
    The frequency is expressed in units of radians/sample.  pmtm uses a 
    default FFT length equal to the greater of 256 and the next power of
    2 greater than the length of X.  The FFT length determines the length
    of Pxx.
...

如果您追求性能，您可能希望使用 2D fft (fft2) 获得信号的频谱表示，然后自行计算功率分布。这将允许您在没有任何 for 循环的情况下处理 2D 数组，但这肯定意味着您需要更多编码:-(

【讨论】：