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

Ver.1 から以下の点について改良した．

<改良点>

• データの先頭と終端まで，フレームの重みを均一化．
• スペクトルとスペクトル平均のプロットに，縦横軸の表記を追加．
• MATLABのエディタで開いたとき，階層化されるようコメントを整理．

<ソースコード>

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

[wavdat, Fs]=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) );

fftf=Fs*( (1:NFFT)-1) / NFFT;

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


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


%% periodgram
tfcomp=zeros(frameN, NFFT);

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

figure;
hoge=pcolor(frameNsec,fftf(1:(NFFT/2+1)),20*log10(abs(tfcomp(:,(1:(NFFT/2+1)))')) );
set(hoge,'EdgeColor','none');
xlabel('Time [sec]');
ylabel('Frequency [Hz]');


%% 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(fftf(1:(NFFT/2+1)),average_spectrum(1:(NFFT/2+1)) - maxdB + maxdBoffset);
xlim([1 fftf(NFFT/2+1)]);
xlabel('Frequency [Hz]');
ylabel('Power [dB]');