Scilabで大容量のCSV(テキスト)ファイルを読み込む



Scilab 5.4.0から、新しい関数が追加されて新しい記事を書きました。下記の内容は古いです。






Scilabを利用される人の中には、データの分析に使う人が多いのではないのでしょうか。
私もその一人で、生体データや振動信号などのデータをScilabに取り込んで、信号処理して表示させるということを良く行います。

データはCSVなどテキストファイルで入力することが多いです。
Matlabの場合は、load関数で大体のファイルは読み込めてしまえますが、Scilabの場合は関数を使い分ける必要があります。

テキストファイルを読み込む関数としては、

・read関数
・fscanMat関数
・read_csv関数
・mfscanf関数

があります。
read関数は最も高速ですが、数字以外のデータが含まれるとエラーがでます。
fscanfMat関数は速度と実用のバランスが取れた関数で、数字以外のデータは自動的にスキップしてくれます。
read_csv関数は、文字列としてcsvを読み込む関数で、低速です。
mfscanf関数は低レベルの関数ですが、複雑なフォーマットを持つファイルを読み込むのに使えます。

各関数の違いを確認するために、実行速度を比較してみます。

a=rand(10000,2);
deletefile('test.csv');
write('test.csv',a);

timer();
b1=read('test.csv',-1,2);
t=timer();
printf("read:%f\n",t);

timer();
b2=fscanfMat('test.csv');
t=timer();
printf("fscanfMat:%f\n",t);

timer();
b3=read_csv('test.csv');
t=timer();
printf("read_csv:%f\n",t);

timer();
fp=mopen('test.csv');
b4=mfscanf(-1,fp,'%f %f');
mclose(fp);
t=timer();
printf("mfscanf:%f\n",t);




実行結果

read:0.062400
fscanfMat:0.280802
read_csv:61.635995
mfscanf:0.312002



read関数は高速なことがわかります。
fscanfMat関数もそこそこ高速です。mfscanf関数はfscanfMat関数とほぼ同等です。
read_csv関数はかなり低速で、大きなファイルを読み込むのには向きません。

Scilab 5.4.0から、新しい関数が追加されて新しい記事を書きました。そちらも参考にしてください。




スポンサーサイト

Scilab 5.2.2 リリース

Scilab 5.2.2がリリースされました。

高速化とバグ修正されたようです。

Scilabで時間周波数解析(スペクトログラム)2

前回、Scilabの時間周波数解析をするプログラムがいまいちだとかきましたが、



Scilabで時間周波数解析するプログラムを作成してみました。



stft.sci

使い方

s:解析対象信号

nfft:FFT長さ

Fs:サンプリング周波数

width:窓幅

noverlap:窓をシフトさせるサンプル数

logflag:ログ表示するかどうか。logflag="l"でログ表示

threshold_min,threshold_max:閾値



function smatrix=stft(s,nfft,Fs,width,noverlap,logflag,threshold_min,threshold_max)

[lhs,rhs]=argn(0);

if ~exists("logflag") then logflag="n"; end;
if rhs < 5 then noverlap=nfft/2; end;
if rhs < 4 then width=nfft;end;
wind=window('hn',width);

nlen=max(size(s));

smatrix=zeros((nlen-width+1)/noverlap,nfft);
for k=1:size(smatrix,1)
s_ind=noverlap*(k-1)+1:noverlap*(k-1)+1+width-1;
s_p=s(s_ind).*wind;
s_p=[s_p(:)' zeros(1,nfft-width)];
smatrix(k,:)=fft(s_p);
end
smatrix=smatrix(:,1:$/2);
f=gcf();
f.color_map=jetcolormap(256);
drawlater();

if exists("threshold_min") then
ind=find(abs(smatrix)<threshold_min);
smatrix(ind)=threshold_min;
end;
if exists("threshold_max") then
ind=find(abs(smatrix)>threshold_max);
smatrix(ind)=threshold_max;
end;

if strcmp(logflag,"n")==0 then

grayplot(((0:size(smatrix,1)-1)*(noverlap/Fs))+(width-1)/2/Fs,(0:size(smatrix,2)-1)/size(smatrix,2)*Fs/2,abs(smatrix));
else
grayplot(((0:size(smatrix,1)-1)*(noverlap/Fs))+(width-1)/2/Fs,(0:size(smatrix,2)-1)/size(smatrix,2)*Fs/2,log10(abs(smatrix)));
end


xlabel('Time[s]');
ylabel('Frequency[Hz]');
drawnow();
endfunction




テストプログラムは以下になります。





exec stft.sci;

fname=SCI+"/modules/sound/demos/"+"chimes.wav";
[y,Fs,bits]=wavread(fname);

data=y(1,1:size(y,2));
t=[0:size(data,2)-1]/Fs;

nfft=512;
width=512;
noverlap=64;

clf;
subplot(211);
dm=max(abs(data));
a=gca();a.data_bounds=[0,-dm;t($),dm];
plot(t,data);
subplot(212);
ret=stft(data,nfft,Fs,width,noverlap,logflag="l",threshold_min=0.0000025);
a=gca();a.data_bounds=[0,0;t($),Fs/2];





test_stft.jpg



良さそうですね。








続きを読む

Scilabで時間周波数解析(スペクトログラム)

Scilabには時間周波数解析をする関数が準備されています。
mapsound という関数です。

試してみましょう。

//テスト用信号をを生成
Fs=22050; //サンプリング周波数
Tl=1; //信号の長さ[s]
t=soundsec(Tl,Fs);

f1=3000; //信号に含まれる最低周波数
s=(sin(2*%pi*t.*(f1+t*f1))); //チャープ信号

scf(); //グラフィックウインドウを開く
mapsound(s,0.02,200,10000,1,Fs); //スペクトログラムをプロット
f=gcf();
f.color_map=jetcolormap(256); //色を変更


mapsound.jpg
表示はできましたが、あまり具合は良くないですね。
オーバーラップさせられない、窓関数かけてないなど、基本的な機能が欠けています。
あくまで簡易的なものと考えたほうが良さそうです。

Scilabで時間周波数解析する関数を自作してみましたので、そちらも見ていただけると幸いです。



scilabでフーリエ解析(簡単版)

前回の記事では、自分でフーリエ変換してプロットする関数を紹介しましたが、
Scilabには最初からフーリエ変換してプロットしてくれる関数が存在します。

analyze という関数がそれです。

試してみましょう。


t=soundsec(1);

freqlist=[100 200 300];
Fs=1/(t(2)-t(1));

sig=sin(2*%pi*t*freqlist(1))+0.7*sin(2*%pi*t*freqlist(2))+0.8*sin(2*%pi*t*freqlist(3));

analyze(sig,10,1000,Fs)


analyze
良さそうです。
縦軸をログにしてみます。

a=gca();
a.log_flags="nln";


analyze2.jpg
良い感じです。

*この記事は、以前のブログの記事scilabでお手軽に周波数解析と時間周波数解析を再編集して再掲しました



scilabでフーリエ解析

フーリエ解析をScilabでしてみます。
fftplot


まずは、フーリエ変換してプロットする関数をつくります。

fftplot.sci
unction [result,freqlist]=fftplot(data,Fs)

fftdata=fft(data);

freqlist=linspace(0,1,length(data))*(Fs-Fs/length(data));

plot2d(freqlist,abs(fftdata),logflag="nl");

a=get("current_axes");
a.data_bounds=[0,min(abs(fftdata)); Fs/2,max(abs(fftdata))];

xlabel('Frequency[Hz]');

result=fftdata;
endfunction



テストプログラム

exec("fftplot.sci");


t=linspace(0,1,1000+1);

freqlist=[100 200 300];
Fs=1000;

sig=sin(2*%pi*t*freqlist(1))+0.7*sin(2*%pi*t*freqlist(2))+0.8*sin(2*%pi*t*freqlist(3));


scf;
subplot(211);
plot2d(t,sig);
ax=get("current_axes");
m=max(abs(sig));
ax.data_bounds=[0 -m; 0.2 m];
subplot(212);
fftplot(sig,Fs);



参考にしたページ
非公式広島大学病院脳磁図室ホームページSCILAB篇
Scilab日本語ドキュメント
Scilab つかいませんか





*この記事は、以前のブログの記事scilabでフーリエ解析を再編集して再掲しました
全記事表示リンク

全ての記事を表示する

最新記事
最新コメント
最新トラックバック
カテゴリ
検索フォーム
RSSリンクの表示
リンク