音の周波数特性を描く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オクターブ幅になっていたのを修正しました。