HHT中toimage函数out of memory问题的解决技巧
先看看函数toimage函数:
[im,tt,ff] = toimage(A,f,t,splx,sply)
Transforms a spectrum made of 1D functions (e.g., output of "hhspectrum") in an 2D image,将1维函数构成的谱变换成2维图像。
A,f,t就是函数"hhspectrum"的输出,output of "hhspectrum",不多说。
splx:number of columns of the output im (time resolution).If different from length(tt), works only for uniform sampling.将时间平分的份数,一般为length(t)。
sply:number of rows of the output E (frequency resolution).
sply表示输出结果E的行数了,即将频率需要平分的份数。
正常频率范围是0~fs/2,如果sply能达到fs/2,那么频率就能分辨出1Hz的值了。
现实中sply选择过大,往往会出现out of memory问题,主要问题就是存储im的矩阵超出matlab内存范围,其实im是一个稀疏矩阵,toimage中计算im是通过下面语句实现:
im = accumarray([indf(:),indt(:)],A(:),[sply,splx]);
如果将上述语句改成:
im = accumarray([indf(:),indt(:)],A(:),[sply,splx],[],[],true);
它表示以稀疏矩阵的形式存储im,因为accumarray函数的完整形式是:
A = accumarray(SUBS,VAL,SZ,FUN,FILLVAL,ISSPARSE)
其最后一个参数就表示可以用稀疏矩阵的形式表示。这样就节省很多存储空间,out of memory问题就能解决。稀疏矩阵im的图标与正常变量的图标不一样,不用管它,以正常矩阵运算即可。
附件是坛友的一份9000*2文件(附件可以点击“阅读原文”登录论坛下载),选用第二列数据,采样点数9000,采样频率500000Hz,我的电脑内存是2G,如果采用原先用法,sply设为10000即报out of memory问题,如果采用改进用法sply设为200000也可正常运行,但设为250000就出现Maximum variable size allowed by the program is exceeded问题。
请大家试试。下面这个仿真信号,将sply设为fs/2也不会出现上述问题,如果不改的话,就出现out of memory。
clc;clear;
fs=102400;N=10240;
n=1:N;t=n/fs;
x1=(1+0.2*sin(2*pi*7.5*t)).*(cos(2*pi*30*t+0.5*sin(2*pi*15*t)));
x2=sin(2*pi*120*t);
x=x1+x2;
imf=emd(x);
[A,fa,tt]=hhspectrum(imf);
[E,tt1,Cenf]=toimage(A,fa,tt,length(tt),fs/2);
nonzero=nnz(E);
%求E中非零项个数。
如果你的电脑内存是4G甚至8G,16G,可能就不会出现上述情况。
声明:本文由声振论坛会员shuihai707原创,转载请注明:来自@声振之家微信公众号。