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

shingoushori's dialy

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

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