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



良さそうですね。








2011.12.4 コピペミスを修正しました
スポンサーサイト

コメントの投稿

非公開コメント

error

ブログの投稿ありがとうございます。今参考にさせて頂いています。
ですが、少しエラーになったところがありますけれども、ここです:
ind=find(abs(smatrix) smatrix(ind)=threshold_min;
どうも、括弧が足りないので、追加しました。
でも、何故か括弧を追加しても、エラーが消えないのです。解決方法は思いつきませんか?
急に質問して、すみません。

ありがとうございます

コピペしたときに、FC2ブログの仕様で、文字化けしていたようです。

修正しました。ごしてきありがとうございます。
全記事表示リンク

全ての記事を表示する

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