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

20 views (last 30 days)
Michael
Michael on 14 Oct 2013
Commented: Michael on 18 Oct 2013
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
  1 Comment
Michael
Michael on 14 Oct 2013
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).

Sign in to comment.

Answers (1)

Wayne King
Wayne King on 14 Oct 2013
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
  2 Comments
Michael
Michael on 18 Oct 2013
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!
take some reading through but the general conclusion seems to be: ' we don't know '...

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by