FFT and Parseval's Theorem

161 views (last 30 days)
Chen Firestein
Chen Firestein on 22 Feb 2021
Edited: David Goodmanson on 5 Mar 2021
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.
  1 Comment
Bjorn Gustavsson
Bjorn Gustavsson on 22 Feb 2021
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...

Sign in to comment.

Answers (1)

David Goodmanson
David Goodmanson on 22 Feb 2021
Edited: David Goodmanson on 5 Mar 2021
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.
  2 Comments
Chen Firestein
Chen Firestein on 5 Mar 2021
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
David Goodmanson on 5 Mar 2021
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.

Sign in to comment.

Tags

Community Treasure Hunt

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

Start Hunting!