周波数応答 / Frequency Response

function [f, MX] = magnitude(x, Fs, octdiv, doMean)
%[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)
%
%   [f, MX] = magnitude(x, Fs, N) can be used to get 1/N-octave band
%   response.
%
%   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
%   2009-09-16 


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

%cfreqs = [12.5, 16, 20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 800 1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500 16000 20000 25000 31500 40000 50000 63000 80000 100000 125000 160000 200000 250000 315000 400000 500000];
if nargin<3
  octdiv = 3;  % 1/3-octave analysis (default)
end
if nargin<4
  doMean = true;   % averaging of magnitude response
end
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) = mean(10*log10(MXtmp(tL(k):tH(k))));
end

% if you need averaging of magnitude response
mxLog = 10*log10(MXtmp);
if doMean
  fnL = 500;
  fnH = 2000;
  [v,indL] = min(abs(ftmp - fnL));
  [v,indH] = min(abs(ftmp - fnH));
  mxMean = mean(mxLog(indL:indH-1));
else
  mxMean = 0;
end

% plot the result
if nargout == 0
  h = semilogx(ftmp, mxLog-mxMean);
  set(h, 'Color', [.6 .6 .6]);
  hold on;
  tt = logspace(log10(min(cfreqs)), log10(max(cfreqs)),200);
  h = semilogx(cfreqs, mxfit-mxMean, 'b.--');
  hold off;
  
  %xticks = [10 20 50 100 200 500 1000 2000 5000 10000 20000 50000 100000 200000 500000 1000000 2000000 5000000];
  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('Power (dB)');
end

f = ftmp;
MX = MXtmp;
if nargin==3
  f = cfreqs';
  MX = (mxfit-mxMean)';
end