音の周波数特性を描くMatlabプログラム

音の周波数特性を描くMatlabプログラムを改善したので、メモ的に置いておきます。今回の目玉は「1/3オクターブごとの平均値を描くようにした」です。500-2000Hz帯の平均値が0dBになるようにスケーリングしています。

function [f, MX] = magnitude(x, Fs)
%[f, MX] = magnitude(x, Fs)
%   Calculate and correctly scale the output of the FFT function to obtain
%   a meaningful power versus frequency plot.  This function plots the
%   result if no output variables are specified.
%
%   Input:
%      x - signal
%      fs - sampling frequency of x in Hz
%   Output:
%      f - Frequency (Hz)
%      MX - Power (simple mean square, not power spectral density)
%
%   Refer to Spectral Analysis section of the Signal Processing Toolbox
%   manual for information on power spectrum and power spectral density.
%
%   2002-12-28 by Marui, Atsushi
%              based on Matlab Tech Note 1702
%             (http://www.mathworks.com/support/tech-notes/1700/1702.shtml)
%   2004-09-20 correspond with the updated Tech Note 1702
%   2005-02-13 added plotting routine
%   2005-08-30 improved plot routine to adjust xlim based on fs
%   2009-07-13 added 1/3-octave line

% Next highest power of 2 greater than or equal to length(x):
NFFT = 2 ^ nextpow2(length(x));

% Take fft, padding with zeros, length(FFTX)==NFFT
FFTX = fft(x, NFFT);
NumUniquePts = ceil((NFFT+1)/2);

% fft is symmetric, throw away second half
FFTX = FFTX(1:NumUniquePts);
MXtmp = abs(FFTX);   % Take magnitude of X

% Scale the FFT so that it is not a function of the length of x.
MXtmp = MXtmp / length(x);
MXtmp = MXtmp .^ 2; % scale properly
MXtmp = MXtmp * 2;  % Account for throwing out second half of FFTX above
MXtmp(1) = MXtmp(1) / 2;                            % Account for DC uniqueness
if ~rem(NFFT,2)
  MXtmp(length(MXtmp)) = MXtmp(length(MXtmp)) / 2;  % Account for Nyquist uniqueness
end
ftmp = (0:NumUniquePts-1) * Fs / NFFT;

% plot the result
if nargout == 0
      octdiv = 3;  % 1/3-octave analysis
      cfreqs = 1000 * 2 .^ ((-20*(octdiv/3):50*(octdiv/3)) / octdiv);
      cfreqs = cfreqs(cfreqs<=Fs/2);
      fcL = cfreqs / 2^(1/(octdiv*2));
      fcH = cfreqs * 2^(1/(octdiv*2));
      tL = zeros(size(cfreqs));
      tH = zeros(size(cfreqs));
      for k=1:length(cfreqs)
        [v,indL] = min(abs(ftmp - fcL(k)));
        [v,indH] = min(abs(ftmp - fcH(k)));
        tL(k) = indL;
        tH(k) = indH-1;
      end
      mxfit = zeros(size(cfreqs));
      for k=1:length(cfreqs);
        mxfit(k) = geomean(MXtmp(tL(k):tH(k)));
      end
      
      fnL = 500;
      fnH = 2000;
      [v,indL] = min(abs(ftmp - fnL));
      [v,indH] = min(abs(ftmp - fnH));
      mx = MXtmp ./ geomean(abs(MXtmp(indL:indH-1)));
      
      [v,indL] = min(abs(cfreqs - fnL));
      [v,indH] = min(abs(cfreqs - fnH));
      mxfit = mxfit ./ geomean(abs(mxfit(indL:indH-1)));
      
      h = semilogx(ftmp, 10*log10(mx));
      set(h, 'Color', [.5 .5 .5]);
      hold on;
      semilogx(cfreqs, 10*log10(mxfit), 'bo');
      semilogx(cfreqs, 10*log10(mxfit), 'b--');
      hold off;
      
      xticks = [16 31.5 63 125 250 500 1000 2000 4000 8000 16000 32000 64000 128000 256000 512000 1024000 2048000];
      k = sum(xticks <= Fs/2);
      xticks = xticks(1:k+1);
      %set(gca, 'XLim', [min(xticks) max(xticks)]);
      set(gca, 'XLim', [20 20000]);
      %set(gca, 'YLim', [-65 15]);
      set(gca, 'XTick', xticks);
      set(gca, 'XTickLabel', xticks);
      grid on;
      set(gca, 'XMinorGrid', 'off');
      set(gca, 'XMinorTick', 'off');
  title('Power Spectrum');
  xlabel('Frequency (Hz)');
  ylabel('Level (dBr)');
else
  f = ftmp;
  MX = MXtmp;
end

【追記 2009-08-06】1/3オクターブ幅でスムージングをするはずが、4/3オクターブ幅になっていたのを修正しました。