99 views (last 30 days)

Show older comments

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
on 22 Feb 2021

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.

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.

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

Start Hunting!