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