読者です 読者をやめる 読者になる 読者になる

shingoushori's dialy

音信号処理を専ら研究していた元博士後期課程の学生によるメモ

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

MATLAB / Octave 信号処理 / Signal Processing

(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);