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

3 次查看（过去 30 天）
Michael
Michael 2013-10-14
评论： Michael 2013-10-18
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 个评论
Michael
Michael 2013-10-14
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).

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Wayne King
Wayne King 2013-10-14
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 个评论
Michael
Michael 2013-10-18
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
Michael 2013-10-18
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 '...

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Estimation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by