关于HHT中toimage函数经常出现out of memory问题的解决技巧

shuihai707

先看看函数[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，可能就不会出现上述情况。

cwb

好厉害啊，这问题困扰我很久了，为此还专门买了2G内存条装上

cwb

搜了一下，觉得 accumarray函数还真是复杂

轻安自在

{:{39}:}{:{39}:}

binglanyouyou

楼主 在disp_hhs出现的这个问题碰到过吗

binglanyouyou

用楼主这个方法改过之后，在函数disp_hhs()时出现溢出内存错误，请问楼主怎么解决啊，非常感谢啊

shuihai707

binglanyouyou 发表于 2014-10-15 21:03
用楼主这个方法改过之后，在函数disp_hhs()时出现溢出内存错误，请问楼主怎么解决啊，非常感谢啊

把数据情况说明清楚，把程序贴一下，这样分析不出来。

binglanyouyou

shuihai707 发表于 2014-10-16 18:21
把数据情况说明清楚，把程序贴一下，这样分析不出来。

clear
clc
close all hidden
sf=12800; %采样频率
%N=40000;
A=xlsread('Book2.xlsx'); %读取excel文件  s');
fid = A(1:10000);
N=length(fid);
t=1/sf:1/sf:N/sf;
figure(1)
plot(t,fid);
%FFT变换
X=abs(fft(fid));
figure(2);
plot((1:N/2)*sf/N,X(1:N/2)*2/N);
%EMD变换
IMF = emd(fid);
[m,n] = size(IMF);
figure(3)
for i = 1:(m-6)
    subplot(m-6,1,i);
    plot(t,IMF(i,1:end));
end
%Hilbert变换
[Aa,fa,tt]=hhspectrum(IMF);
[im,tt1]=toimage(Aa,fa,tt,length(tt));
disp_hhs(im);
for k=1:size(im,1)
    bjp(k)=sum(im(k,:))*(1/sf)*(sf/N);
end
figure(4);
f2=(0:N-3)/N*(sf/2);
plot(f2,bjp);  % 作边际谱图   进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率，所以此时的频率轴已经是中心频率
xlabel('频率 f/Hz');
ylabel('幅值');
title('边际谱');我的程序是从论坛下载的，已附上，数据是别人测量的振动加速度数据，分析时disp函数总是出现溢出错误，尝试着分析前一万个点，结果仍然溢出，无法求出Hilbert谱以及边际谱

shuihai707

binglanyouyou 发表于 2014-10-18 14:40
clear
clc
close all hidden

disp_hhs函数的输入项有一个是inf (optional): -dynamic range in dB (wrt max)，图像中右侧的能量条以dB的形式显示，实际上我们完全可以用其他形式取代，直接调用imagesc(tt/sf,0:0.5*sf,im)这个函数也能实现Hilbert谱画图，这样log2内存出错就不会出现，但是由于横、纵坐标范围都很大，也比较占内存。disp_hhs函数本质上也是调用这个函数，它将im的能量通过语句M=max(max(im));im = 10*log10(im/M);转化成dB的形式了。至于为什么log2函数对稀疏矩阵求对数运算会出现内存溢出问题，我也不是很清楚。

binglanyouyou

shuihai707 发表于 2014-10-18 22:36
disp_hhs函数的输入项有一个是inf (optional): -dynamic range in dB (wrt max)，图像中右侧的能量条以dB的 ...

imagesc(tt/sf,0:0.5*sf,im)，您好 这句程序是您实验过可以运行的吗 为什么我的还是不可以

shuihai707

binglanyouyou 发表于 2014-10-20 10:48
imagesc(tt/sf,0:0.5*sf,im)，您好 这句程序是您实验过可以运行的吗 为什么我的还是不可以

如果还是内存不足，可能就是矩阵太大的原因，我昨天试了一下，10000个数据，我2G的内存，程序不报错，但图像出不来，不知何故，程序运行中内存突增。

binglanyouyou

shuihai707 发表于 2014-10-20 15:36
如果还是内存不足，可能就是矩阵太大的原因，我昨天试了一下，10000个数据，我2G的内存，程序不报错，但 ...

谢谢您这么费心，真是很头疼，您能把您在处理实际数据的程序发给我一份吗 我的是不是程序的问题呢

shuihai707

binglanyouyou 发表于 2014-10-21 10:04
谢谢您这么费心，真是很头疼，您能把您在处理实际数据的程序发给我一份吗 我的是不是程序的问题呢

你的程序没有问题，我就是用你的程序仿真的，我也发帖求助了，有解决方法告诉你

shuihai707

binglanyouyou 发表于 2014-10-21 10:04
谢谢您这么费心，真是很头疼，您能把您在处理实际数据的程序发给我一份吗 我的是不是程序的问题呢

sf=12800; %采样频率
fid = data(1:10000);
fid=fid-mean(fid);
N=length(fid);
t=1/sf:1/sf:N/sf;
IMF = emd(fid);
[Aa,fa,tt]=hhspectrum(IMF);
[im,tt1,Cenf]=toimage(Aa,fa,tt,length(tt),sf/2);
for k=1:size(im,1)
    bjp(k)=sum(im(k,:))*1/sf;
end
figure;
plot(Cenf(1,:)*sf,bjp);xlim([0 200]);
xlabel('频率 f/Hz');
ylabel('幅值');
title('边际谱');