MATLAB: Energy of a signal in t and f domain

digital signal processingenergyfftsignalspectrum

The energy of a signal is expected to be the same in t and f domain.

n  = 1e4;
dx = 0.25;
x  = rand(n,1) -0.5;
ex = sum(x.^2) *dx; % energy in t domain
y  = fft(x);
fs = 1/dx;
df = fs/n;
ya = abs(y);
ey = sum(ya.^2) *df; % energy in f domain

but from the code, ey/ex=16, exactly the squared fs.

what's the problem?

Best Answer

I found the solution myself.

Normalizing spectral amplitude by fs will work.

But I don't know why to do this.

Related Solutions

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

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.

MATLAB: Code for ploting fourier transform spectrum

spectrogram()

Best Answer

Related Solutions

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

MATLAB: Code for ploting fourier transform spectrum

Related Question