MATLAB: Physically correct Normalization of fft + Implementing Parseval’s Theorem

fftMATLABparseval's theoremscaling fft

Hi. I've read endless amounts of similar stated questions on this forum, but I can't seem to find the answer to my specific question.
I would like to tranform a timesignal recorded by an accelerometer into the frequency domain. The accellerometer measured the acceleration of a hammer during impact with a floor. The Idea is to use the acceleration to calculate the force of the impact and then use the results to perform a FEM modelling of a floor. To save computation time in the FEM-model I wish to do the calculations in the frequency domain, since I'm not interested in the transient. Now, my question is as follows:
How should I scale the frequency amplitudes resulting from the fft, to get physically correct amplitudes in the frequency domain? After searching the forum and elsewhere for days I've come up with the following script:
%FFT 
[x,Fs] = audioread('acceleration-file.WAV'); %x = Signal, Fs = samplingfrequency
T=1/Fs;             %Sample period
L=length(x);        %Signal length
df=2*Fs/L;
% FFT
NFFT=2^nextpow2(L);
Y  = fft(x,NFFT)*T;                 %FFT of signal--> Will the multiplication with T result in a physically correct amplitude?
Ya = abs(Y);                        % Normalization of amplitude
f = Fs/2*linspace(0,1,L/2+1);       %Frequency vector, only positive frequencies
YA= 2*Ya(1:L/2+1);                  %Frequency amplitudes, times two for one-sided spectrum
YA(1)=YA(1)/2;                      %Correction for the fact that DC doesn't occur twice                      
%Plot one-sided amplitude spectrum
figure
plot(f,YA)
title('Onesided amplitude spectrum')
xlabel('Frequency [Hz]')
ylabel('|y(t)|')
grid on
%Checking Parseval's Theorem

energy_x=sum(x.*conj(x)*T);
energy_YA=sum(YA.*conj(YA)*df);
diff=energy_x-energy_YA            %Doesn't return 0.
Is this normalization correct? Also, is the implementation of Parseval's Theorem correct? To be clear, it's the following lines I'm really curious about(since parseval's theorem doesn't seem to hold):
%Normalizing amplitudes  
Y  = fft(x,NFFT)*T;
%Checking Parseval's Theorem
energy_x=sum(x.*conj(x)*T);
energy_YA=sum(YA.*conj(YA)*df);
diff=energy_x-energy_YA
Any help would be greatly appreciated, as satisfactory explanations are hard to come by out there.

Best Answer

Hello fgb
Here is a revised answer that attempts to better address the issue. I believe there are a couple of things going on here. In your code,
   df = 2*Fs/L
should not have the factor of 2. Also, if the length of the x array is not a power of 2, then the procedure
NFFT = nextpow2(L)
Y ~ fft(x,NFFT)
changes the length of the x array that gets fft'd. So both uses of L afterwards have to be changed, and df has to become
df = Fs/NFFT.
For a transient response the transform and the energies in the time and frequency domains are
Y = fft(x)*dt                % dt is the same as T
ex = sum(conj(x).*x)*dt      % dt is the same as T 
eY = sum(conj(Y).*Y)*df
The first two you have. The frequency domain calculation eY agrees with ex, and it uses the entire frequency array. All frequency components have identical energy content. Each energy is the abs value squared of the amplitude. Life is good. The issue is with using positive frequencies only.
When you fft a real quantity and double the positive positive frequency components Y(f), you are going to a representation (just at one frequency for simplicity)
x = 2 |Y(f)| cos(2pi f t + phase angle)
But the energy content is different, because the average value of cos(x)^2 is 1/2. So if somebody has handed you a positive-frequency-absolute-value spectrum like the one that you have after this step
YA(1)=YA(1)/2;  %Correction for the fact that DC doesn't occur twice
then you also have to make sure that point N/2+1 has been included (as you did), and that that point is not doubled. So there needs to have been a step
YA(end)=YA(end)/2;  %Correction because nyquist point doesn't occur twice
Then the expression
eY_new = (YA(1)^2 + (1/2)*sum(YA(2:end-1).^2) + YA(end)^2)*df;
works. But it's not straightforward.
Stuff like this is why I try to stay away from positive-frequencies-only as much as possible for real calculations. It's not bad for illustrative plotting purposes.
Related Question