FFT amplitude scaling - sample frequency/sampling period, number of input samples, transform length, or none of the above?

I seem to be reading conflicting information from different sources around Matlab Central on this. How should the amplitude scaling of a signal analysed using FFT be carried out?

I have analysed acceleration time-history data using 1/3-octave filtering and FFT, and (assuming the 1/3-octave filtering is correct, performed using the '<http://www.mathworks.com/matlabcentral/fileexchange/69-octave octave>' toolbox) the amplitude scaling is, unsurprisingly, different. From what I have read I need to either divide the FFT output by the number of signal samples (the length of the input vector) or the transform length (the length of the output vector). However, when I do this, the amplitude scaling is different by a factor of around 10 for each (since in my case the length of the input and output vectors differ only by around 30% there is little relative difference between the two). The following are image files of the output graphs with the y-axis amplitude scales shown.

My code looks like this:

      Fref = [2.5, 3.1 5 4 5, 6.3 8 10, 12.5 16 20, 25 31.5 40, 50 63 80];  %    Preferred labeling freq. from BS EN ISO 266:1997
      n = -26: 1: -11;
      Fc = 1000*(10.^(n/10));        %    Exact center freq. from BS EN ISO 266:1997  

      %   Fast Fourier Transform

      m = length(samples);          % Window length
      nfft = pow2(nextpow2(m));  % FFT transform length
      P1fft = fft(P1,nfft);     % FFT
      P1fft = P1fft(1:nfft/2)/length(samples);  % Discard unnecessary part of spectrum and scale amplitude

      P2fft = fft(P2,nfft);     % FFT
      P2fft = P2fft(1:nfft/2)/length(samples);  % Discard unnecessary part of spectrum and scale amplitude

      f = (0:nfft/2-1)*(fs/nfft);     % Frequency range

      figure
      semilogx(f,abs(P2fft),f,abs(P1fft))
      grid
      xlim([min(Fref) max(Fref)])
      % ylim([0 0.014])
      set(gca,'XTick',Fref)
      xlabel('Frequency (Hz)')

      %   Preallocation for filtered data matrices

      P1filt = zeros(length(samples),length(Fref));
      P2filt = zeros(length(samples),length(Fref));

      %   Filter design

      N = 3;  %   order of filter

      figure
      for i = length(Fc):-1:1 %   Backwards indexing equivalent to preallocation
        [B,A] = oct3dsgn(Fc(i),fs,N);
        y(:,i) = filter(B,A,P1);
        P1filt(:,i) = y(:,i);
        P1rms(:,i) = sqrt(mean(P1filt(:,i).^2));

              [G,F] = oct3spec(B,A,fs,Fc(i),'ANSI',N);  %   Display 1/3-octave filter spec
              semilogx(F,G)
              ylim([-100 3])
              hold on

        [B,A] = oct3dsgn(Fc(i),fs,N);
        y(:,i) = filter(B,A,P2); 
        P2filt(:,i) = y(:,i);
        P2rms(:,i) = sqrt(mean(P2filt(:,i).^2));
      end

      figure
      semilogx(Fref,P2rms,Fref,P1rms)
      grid
      xlim([min(Fref) max(Fref)])
      set(gca,'XTick',Fref);    % Label frequency axis.
      xlabel('Frequency (Hz)')

The input signal vectors, 'P1' and 'P2', are acceleration time-history waveforms. The variable 'samples' is also a vector of the acceleration values - the reason for having an independent variable containing the same information as the input vectors is because I'm processing multiple waveforms, that all have the same length as 'samples'.

Please could any helpful contributors advise on correct FFT amplitude scaling, and hopefully put this to bed for other users also?!

Any general advice on optimising or improving my code would also be much appreciated!

Cheers,

Mike

[Michael]

I want to clarify: the factor 10 difference found applies to fft division by either of the transform length or input signal length (i.e. it is not a factor of 10 difference between on the one hand, dividing by the transform length, and on the other, the input signal length).

Answer 1

[Wayne King]

Do you have the Signal Processing Toolbox in 13a or 13b? If so, then use periodogram() with the 'power' option. That will do things correctly for you.

    Fs = 1000;
    t = 0:1/Fs:1-1/Fs;
    x = cos(2*pi*100*t);
    [Pxx,F] = periodogram(x,[],length(x),1000,'power');
    plot(F,Pxx)

Note the power estimate is 1/2, which is exactly what you would expect.

Now see how it stays correct even if I change the FFT length and use a window.

      x = 2.5*cos(2*pi*100*t);
     [Pxx,F] = periodogram(x,flattopwin(length(x)),1024,1000,'power');
      val = max(Pxx)

Note the power estimate is 3.125, which is exactly what you expect (2.5)^2/2

[Michael]

Thanks for your response Wayne. I do have that toolbox, but I'm in 12a and the 'power' option doesn't appear to be available. I'm also not looking for the PSD or power, but am seeking to compare acceleration magnitude spectra directly. Shouldn't the standard fft function be able to do this, with appropriate scaling?

[Michael]

Wow, I've read through some of the other posts on the subject scattered about Matlab Central and this seems to have been a contentious issue for some time!

http://www.mathworks.co.uk/matlabcentral/answers/39344 http://www.mathworks.co.uk/matlabcentral/newsreader/view_thread/80750

take some reading through but the general conclusion seems to be: ' we don't know '...