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