MATLAB / Octaveで，スペクトログラムやスペクトル平均を算出 Ver.1

(2013.08.20 追記) 改良版 Ver.2をつくりました -> Ver.2

spectrogram関数を用いてしまえば，一発でプロットできてしまいます．

短時間フーリエ変換を使用したスペクトログラム - MATLAB spectrogram

しかし，窓関数やゼロ詰めといった諸元を調整するには，手間がかかります．
ならばいっそ，FFT以外は自分で書いてしまおうという試み．

<ソースコード>

	諸元
FFT点数	4096点
解析窓	2048点 ハニング窓
フレームシフト	512点
正規化	入力信号を最大振幅値の0 Hzとしたとき，0 Hzが96 dB

wavdat=wavread('hoge.wav');

% setup
NFFT=4096;
windowsize=2048;
framelate=512;

h=zeros(1,NFFT);
h(1:windowsize)=0.5*(1 - cos(2*pi* ( (1:windowsize)-1)/windowsize) );

maxdB=20*log10(sum(h));
maxdBoffset=96;


frameN=ceil( (length(wavdat)+framelate+(NFFT-framelate))/framelate );
newlength=frameN*framelate + (NFFT - framelate);
buf_wav=zeros(1,newlength);
buf_wav( (framelate+1):(framelate+length(wavdat)))=wavdat;


% periodgram
tfcomp=zeros(frameN, NFFT);

for i=1:frameN
    tfcomp(i,:)=fft( h .* buf_wav( (1:NFFT)+(i-1)*framelate) );
end

figure;
hoge=pcolor( 20*log10(abs(tfcomp(:,(1:(NFFT/2+1)))')) );
set(hoge,'EdgeColor','none');


% average spectrum
average_spectrum=zeros(1,NFFT);

for k=1:NFFT
    average_spectrum(1,k)=10*log10( sum( abs(tfcomp(:,k)).^2 / frameN ) );
end

figure;
plot(average_spectrum - maxdB + maxdBoffset);