FFT and Parseval's Theorem

When I try to implement Parseval's theorem, I find that results are incorrect in the frequency domain. Here, my code is attached .

t=linspace(0,10,1000);
Ts=t(2)-t(1);
AA=cos(2*pi*t);

TD_Result=trapz(t,abs(AA).^2);


frequencyAxis = ((0:length(t)-1) -ceil((length(t)-1)/2))/length(t)/(t(2)-t(1));
Fs=frequencyAxis(2)-frequencyAxis(1);
AA_S=fftshift(fft(AA));
AA_S=AA_S/(length(t));
Freq_Result=trapz(frequencyAxis,abs(AA_S).^2);

Is the normalization correct? Looking on the spectrum provides amplitude of 0.5 which is correct.

[Bjorn Gustavsson]

You'll have to check the documentation of fft to figure that out. As far as I understand there are at least 3 different conventions for the normalization of FFTs, and two for the sign of the wt-exponent in the FFT. These conventions lead to different normalization-factors for Parseval...

Answer 1

[David Goodmanson]

Hi Chen,

Here are a couple of ways. Since the time function fills up the entire time window (as opposed to being some kind of pulse that dies off before reaching each end), it is a bit easier to get a close equality by using the rectangular rule approximation for each integral, [sum of points] x [width of interval]. The trapezoidal rule is not quite right for a periodic window.

You are using Fs for the frequency array step width, which conflicts a bit with the usual notation of Fs as the sampling frequency, but I kept Fs as you have it.

By introducing a time array and trapz, you are implicitly declaring the fft to be an approximate representation of a fourier transform integral. Accordingly in the expression for AA_S, there is a factor of Ts for the width of each interval in the riemann sum.

The TD and Freq results below agree.

N = 1000;
t = ((0:N-1)/N)*10;    % time record is not quite to 10 but new time interval is exactly .001
Ts=t(2)-t(1);
AA=cos(2*pi*t);
TD_Result=sum(abs(AA).^2)*Ts        % factor of Ts for width 
Fs = 1/(N*Ts);                      % golden rule for fft
frequencyAxis = ((0:N-1)/N)*Fs;     % not actually used
AA_S=fftshift(fft(AA))*Ts;          % factor of Ts for width
Freq_Result=sum(abs(AA_S).^2)*Fs    % factor of 'Fs' (freq step width) for width

Another way, which dispenses with interval widths and is more signal processing oriented, is

N = 1000;
r = rand(1,N);
E1 = sum(r.^2)
g = fft(r);
E2 = sum(abs(g).^2)/N

When you prove parseval's theorem and plug in ffts, there is a sum over the product of a couple of complex exponentials, and that sum is zero except for one instance where the product of the exponentials is 1. Then the sum over points gives N, which gets compensated for by the 1/N factor on the last llne. Needless to say if you reverse engineeer your way out of the first method above, you can obtain the second method.

[Chen Firestein]

@David Goodmanson

I have another simillar issue: consider the following code for a pulse

Temporal_Pulse_Width=1;
t=linspace(0,20,2000);
dt=t(2)-t(1);
Time_Period=round(Temporal_Pulse_Width/dt);
AA=zeros(length(t),1);
AA(Time_Period/2:3*Time_Period/2)=1;
A=20*log10(abs(fftshift(fft(AA/length(t)))));
frequencyAxis = ((0:length(t)-1) -ceil((length(t)-1)/2))/length(t)/(t(2)-t(1));
plot(frequencyAxis,A,'LineWidth',3);
max(A)

The result (max(A)) is changed when t is changed, i.e. if you take instead

t=linspace(0,10,2000);

The result will be different. How can I do good normalization to my fft which keeps the amplitude in that case?

Thanks,

Chen

[David Goodmanson]

Hi Chen,

Looks like you basically want to do a fourier integral. Try normalizing as in the first method in the answer, and see what happens.