Why isn't the energy preserved in FFT?
Gui Xie
By Parsevel's theorem, the energy of a signal expressed in the time
domain to the energy expressed in the frequency domain should be equal.
However, this cannot be verified by the following Matab code:
F=fft(S); % S is a given signal
E1=sum(S.^2); %the energy in the time domain
E2=sum(abs(F).^2); %the energy in the frequency domain
E2 is much larger than E1. What's wrong with the FFT implementation in
Matlab?
Athanasios Doukas
I think that this happens because FFT is already a sum (see the help)
that takes into account all the other parts of the signal in order to
give you the result. So i do not think that this is a bug or
something like that.
I hope that i helped you.
César Delgado González
I think the time domain energy can't be calculated as a sum of
coefficients. On frecuency domain, this is right cause the FFT give
you the Fourier Expansion coefficents, so energy, on frecuency
domain, could be calculate using Parseval Theorem:
Ef = sum(abs(Xf).^2);
But on time domain you should do the integration of the
abs(signal). The result should be something related to the variance
of the signal on this way:
Et = k*(std(signal).^2);
I don't remember the k value, see for gaussian distribution this
value.
Hope this helps,
Cėsar
Rune Allnor
It's a matter of implementation. From the FFT help-file:
For length N input vector x, the DFT is a length N vector X,
with elements
N
X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
n=1
The inverse DFT (computed by IFFT) is given by
N
x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
k=1
The key is the 1/N factor in the IFFT. If you look up the texts on
the DFT, you will find that both the FFT and the IFFT are scaled by
a factor 1/sqrt(N). I don't know why, but I have seen the same thing
in a number of implementations of the FFT/IFFT; that the FFT is not
scaled at all while the IFFT is scaled by 1/N.
For Parseval's theorem to hold with data, you need to do something
like this:
s=randn(N,1);
S=1/sqrt(N)*fft(s);
es=s'*s;
ES=S'*S;
and the difference es-ES vanishes.
Rune
Derek Goring
If only life were that simple, we'd all be professors of mathematics.
Try this:
S=randn(1000,1);
and run your routine, then calculate
E2/E1
Now try it again with:
S=randn(2000,1);
and calculate
E2/E1
Do you see a pattern emerging?
César Delgado González
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all, clc, close hidden all;
% 1.- CALCULATION
%-----------------
N = 1000; % Number of iterations
for k=1:N
s = randn(k*10,1);
S = fft(s);
Et1(k) = (10*k.*std(s)).^2;
Et2(k) = 10*k*sum(s.^2);
Ef(k) = sum(abs(S).^2);
s=[]; S=[];
end;
% 2.- PLOTTING RESULTS
%---------------------
% 2.1.- Time Variance versus Freq. Sum
figure(1)
plot(Et1,Ef,Et1,fix(mean(Et1)).*Et1./Ef);
title('Energy calculation');
xlabel('Time domain as variance');
ylabel('Frecuency domain');
grid on; zoom on;
% 2.1.- Time sum versus Freq. Sum
figure(2)
plot(Et2,Ef,Et2,fix(mean(Et2)).*Et2./Ef);
title('Energy calculation');
xlabel('Time domain as sum');
ylabel('Frecuency domain');
grid on; zoom on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
So, the k factor is length(signal). I think that this is more
concerned with the variance calculation process more than the base
normalization one.
All bases are usually normalized when are used. Basically two are
the ways to do that:
(1/sqrt(T))*base <------> (1/sqrt(T))*base
(1/T)*base <------------> base
On Fourier, T=2*pi, for that reason, on one way you see
1/(2*pi),and on the other way does not appear this term. You can do
the same using, on both ways, 1/sqrt(2*pi).
Probably the N factor on one way of FFT is in order to avoid the
sqrt operation, just for speed reasons.
So, I think that the length(signal) term is more concerned with
variance calculation proccess than base normalization, that always is
applied (think on rectangular base, all unitary vector has energy =
1).
Hope this helps,
César
AJ "no z" Johnson
the comb-line-like spectrum. Therefore there's infinite energy. What you are
computing when you do an FFT is computing the power (energy per unit time)
spectrum. [At least this is one interpretation.]
(2) The FFT total power should be scaled by 1/N. Then you'll find E1==E2/N.
-aj
Greg Heath
news:cgvir5$4e...@cui1.lmms.lmco.com...
Yes. I find it useful to use the notation PSD for power spectral density.
Then
P1 = sumsqr(f);
F = fft(f);
PSD = abs(F).^2/N;
P2 = sum(PSD);
Hope this helps.
Greg
AJ "no z" Johnson
news:2ph2fkF...@uni-berlin.de...
> "AJ "no z" Johnson" <aj.jo...@lmco.com> wrote in message
> news:cgvir5$4e...@cui1.lmms.lmco.com...
> > "Gui Xie" <g-...@jaist.ac.jp> wrote in message
> > news:cgugn8$km7$1...@jaist-news.jaist.ac.jp...
> > > Hi:
[snip]
> Yes. I find it useful to use the notation PSD for power spectral density.
> Then
>
> P1 = sumsqr(f);
>
> F = fft(f);
> PSD = abs(F).^2/N;
>
> P2 = sum(PSD);
>
> Hope this helps.
>
> Greg
You got me thinking: it is probably computationally cheaper to compute the
squares as follows:
P2 = (1/N)*sum(F.*conj(F));
or
P2 = (1/N)*F*F'; % not sure which (F) the (') goes on
rather than any abs().^2 operation, since that squares a square root.
-Aj
César Delgado González
forever, hence the comb-line-like spectrum. Therefore there's
infinite energy. What you are computing when you do an FFT is
computing the power (energy per unit time) spectrum. [At least this
is one interpretation.
zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
I am not agree with this point. Yes, FFT consider a infinite
periodic signal, that are power signal and not energy signals.
Power signals do not have Fourier Transform, cause Dirichlet
conditions are not achived. But, if the power signal is periodic,
then you can aproximate it with a Fourier expansion, that results on
a sampled spectrum version of one period energy signal.
But this is not the point, the point is that both Fourier
transform and FFT compute J/Hz or J*sg. The power is only achived
integrating a limite bandwidth.
Parseval says that the power of the signal (std^2 over a
bandwidth) is the same as the sum of the Fourier coeficients. So, the
sum over one period of the spectrum is the same as energy over one
period of the periodic time domain signal.
About Power Spectral Density, I don0't like to use this term on
those signals.
All real world signals are stochastical proccess. Only
stacionaries ones has Power Spectral Density, which is the Fourier
transform of the autocorrelation of the signal. This is not the FFT,
is FFT(xcorr(x,x)) on stacionary proccess (every instant of time has
the same stadistical distribution)
But signals are not stacionary, only if we consider small porcion
of signal, where the stadistical parameters haven't got enought time
fot changing.
Then appear many tarnsformation for signal characterization, as
Wavelet, Garbor, etc.
So, the point is:
1.a) PSD = FFT(xcorr(x,x))
1.b) Fourier transform give us the energy distribution over
the frecuency axe. If we integrate a Bw, we get the Power.
This is exactly what we are doing by summing all coeficients
of one period of the spectrum, or calculating the variance of one
period of the time domain signal.
zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
(2) The FFT total power should be scaled by 1/N. Then you'll find
E1==E2/N.
-aj
zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
I think that, probably, this is cause by the integration proccess.
On many time we have to add some multiplication term:
int(exp(-f))= -(1/f)*exp(f) + K
So, probably the length(Signal) factor is concerned with this
proccess.
César
Greg Heath
> "Greg Heath" <he...@alumni.brown.edu> wrote in message
> news:2ph2fkF...@uni-berlin.de...
> > "AJ "no z" Johnson" <aj.jo...@lmco.com> wrote in message
> > news:cgvir5$4e...@cui1.lmms.lmco.com...
> > > "Gui Xie" <g-...@jaist.ac.jp> wrote in message
> > > news:cgugn8$km7$1...@jaist-news.jaist.ac.jp...
> > > > Hi:
> [Parsevel's theorem applies to FFT question]
> [snip]
> > Yes. I find it useful to use the notation PSD for power spectral
density.
> > Then
> >
> > P1 = sumsqr(f);
> >
> > F = fft(f);
> > PSD = abs(F).^2/N;
> >
> > P2 = sum(PSD);
>
> squares as follows:
> P2 = (1/N)*sum(F.*conj(F));
> or
> P2 = (1/N)*F*F'; % not sure which (F) the (') goes on
> rather than any abs().^2 operation, since that squares a square root.
Interesting question. More interesting answer.
I calculated power 7 different ways over M=10 trial runs using N = 2^16
data points. The results varied from trial to trial. Running R13 on WXP, I
got (approx) 1 of 3 answers: Either 0.032, 0.016 or 0.0 seconds.
Code and tabulations for two trials of 10 are given below.
I wouldn't draw any definitive conclusions until more combinations of
M and N were considered.
Calculation Counts(0.032,0.016,0.0)
-------------------------------- ----------------------------
sumsqr(x) (3, 4, 3) & (1, 2, 8)
sum(x.*x) (0, 0, 10) & (0, 1, 9)
sum(abs(X).^2)/N (0, 4, 6) & (1, 8, 1)
X0=abs(X); sum(X0.^2)/N (0, 10, 0) & (0, 9, 1)
X0=abs(X); sum(X0.*X0)/N (0, 5, 5) & (0, 8, 2)
sum(X.*conj(X))/N (0, 6, 4) & (0, 7, 3)
X'*X/N (1, 2, 8) & (2, 5, 3)
clear, close all, clc, format long g
M = 10, N = 2^16
dt = 2*pi/N, t = dt*(0:N-1)';
df = 1/(N*dt), f = df*(0:N-1)';
x = cos(t).*cos(20*t);
%figure, plot(t,x)
X = fft(x);
%figure, plot(f,abs(X).^2/N)
for i = 1:M
tic
P1 = sumsqr(x)
T(i,1) = toc;
tic
P2 = sum(x.*x)
T(i,2) = toc;
tic
P3 = sum(abs(X).^2)/N
T(i,3) = toc;
X0=complex(0,0)*zeros(N,1);
tic
X0 = abs(X);
P4 = sum(X0.^2)/N
T(i,4) = toc;
tic
X0 = abs(X);
P5 = sum(X0.*X0)/N
T(i,5) = toc;
tic
P6 = sum(X.*conj(X))/N
T(i,6) = toc;
tic
P7 = X'*X/N
T(i,7) = toc;
end
sumT = sum(T);
[T;sumT]
Hope this helps,
Greg
AJ "no z" Johnson
[benchmark program]
> Greg
Here is my shot at it, results from 600MHz Windows 2000 machine.
25x improvement, better than I would have guessed.
-Aj
RESULTS:
M = 1000
N = 65536
t= 3.205 X'*X.*(1/N)
t= 5.778 sum(x.*x)
t=77.802 sum(abs(X).^2)/N
t=79.845 X0 = abs(X);P4 = sum(X0.^2)/N;
t=79.885 X0 = abs(X);P5 = sum(X0.*X0)/N
t=19.748 sum(X.*conj(X))/N
t= 3.115 X'*X/N
PROGRAM:
function square_benchmark
% Original by Greg Heath
M = 1000, N = 2^16
dt = 2*pi/N; t = dt*(0:N-1)';
df = 1/(N*dt); f = df*(0:N-1)';
x = cos(t).*cos(20*t);
%figure, plot(t,x)
X = fft(x);
%figure, plot(f,abs(X).^2/N)
pause(0); % Force display update
eq = {};
T = [];
eq{end+1} = 'X''*X.*(1/N)';
tic
for i = 1:M
P1 = X'*X.*(1/N);
%P1 = sumsqr(x) % don't know this function
end
T(end+1) = toc;
eq{end+1} = 'sum(x.*x)';
tic
for i = 1:M
P2 = sum(x.*x);
end
T(end+1) = toc;
eq{end+1} = 'sum(abs(X).^2)/N';
tic
for i = 1:M
P3 = sum(abs(X).^2)/N;
end
T(end+1) = toc;
eq{end+1} = 'X0 = abs(X);P4 = sum(X0.^2)/N;';
tic
for i = 1:M
X0 = abs(X);
P4 = sum(X0.^2)/N;
end
T(end+1) = toc;
eq{end+1} = 'X0 = abs(X);P5 = sum(X0.*X0)/N';
tic
for i = 1:M
X0 = abs(X);
P5 = sum(X0.*X0)/N;
end
T(end+1) = toc;
eq{end+1} = 'sum(X.*conj(X))/N';
tic
for i = 1:M
P6 = sum(X.*conj(X))/N;
end
T(end+1) = toc;
eq{end+1} = 'X''*X/N';
tic
for i = 1:M
P7 = X'*X/N;
end
T(end+1) = toc;
% summarize
for ii=1:length(eq)
fprintf('t=%6.3f %s\n', T(ii), eq{ii});
end
ajab
also a suggestion about this problem:
I've scaled the data by multiplication with frequency matrix, which
is also scaled by sampling rate. But comparing to the equivalent
spatial and frequency operations the magnitudes are not the same!
Below I have coded for 2 case, one calculating the laplacian and the
other is differentiation in x-direction!
% frequency matrix
G=(0:6); H=(-5:-1);
g=repmat(G,12,1); h=repmat(H,12,1);
k=cat(2,g,h);
l=k';
delta=0.1; % sampling interval
fs=2*pi/delta % sampling frequency
% frequency matrix scaled with fs in x- and y-direction)
k=fs.*k;
l=fs.*l;
K=sqrt(k.^2+l.^2); % radial frequency
s=zeros(12); s(2:10,5:7)=1; % signal
S=fft2(s);
lapFreq=real(ifft2(S.*(K.^2))); % laplacian with fft
lap=(del2(s))./fs; % laplacian with 'del2'
dyFreq=real(ifft2(S.*(i*k))); %case2, differentiation with respect to
x by fft
dy=((diff(s'))'./fs); % ... with diff in time domain
subplot(221), imagesc(lapFreq), colorbar, title('lapFreq')
subplot(222), imagesc(lap), colorbar, title('lap')
subplot(223), imagesc(dxFreq), colorbar, title('dxFreq')
subplot(224), imagesc(dx), colorbar, title('dx')
Greg Heath
> "Greg Heath" <he...@alumni.brown.edu> wrote in message
> news:2phvl7F...@uni-berlin.de...
> [benchmark program]
> > Greg
>
> Here is my shot at it, results from 600MHz Windows 2000 machine.
> 25x improvement, better than I would have guessed.
Only if you take the formation of abs(X) into account. Since it is
usually available even when power isn't calculated, I've added
two more calculations which exclude the time for X0=abs(X).
This reduces the ratio by a factor of 8.
> RESULTS:
> M = 1000
> N = 65536
> t= 3.205 X'*X.*(1/N)
> t= 5.778 sum(x.*x)
> t=77.802 sum(abs(X).^2)/N
> t=79.845 X0 = abs(X);P4 = sum(X0.^2)/N;
> t=79.885 X0 = abs(X);P5 = sum(X0.*X0)/N
> t=19.748 sum(X.*conj(X))/N
> t= 3.115 X'*X/N
>> help sumsqr
-----SNIP
SUMSQR(M) returns the sum of the squared elements in M.
(Obviously not the same M as above)
CALCULATION TIME
M=100 M= 1000
-------------------------------- ----------------------------
sumsqr(x) 0.205 1.480
sum(x.*x) 0.189 1.790
sum(abs(X).^2)/N 0.970 9.632
X0=abs(X); sum(X0.^2)/N 0.984 9.992
X0=abs(X); sum(X0.*X0)/N 1.048 10.222
sum(X.*conj(X))/N 0.545 5.922
X'*X/N 0.031 0.434
sum(X0.^2)/N 0.124 1.184
sum(X0.*X0)/N 0.126 1.392
When the creation of abs(X) is taken into account.
0.970/0.031 = 31.29
9.632/0.434 = 22.19
otherwise (since it usually available anyway)
0.124/.031 = 4.00
1.184/0.434 = 2.73
Hope this helps.
Greg
Greg Heath
> "Greg Heath" <he...@alumni.brown.edu> wrote in message
> news:2phvl7F...@uni-berlin.de...
> [benchmark program]
> > Greg
>
> Here is my shot at it, results from 600MHz Windows 2000 machine.
> 25x improvement, better than I would have guessed.
Only if you take the formation of abs(X) into account. Since it is
usually available even when power isn't calculated, I've added
two more calculations which exclude the time for X0=abs(X).
This reduces the ratio by a factor of 8.
> RESULTS:
> M = 1000
> N = 65536
> t= 3.205 X'*X.*(1/N)
> t= 5.778 sum(x.*x)
> t=77.802 sum(abs(X).^2)/N
> t=79.845 X0 = abs(X);P4 = sum(X0.^2)/N;
> t=79.885 X0 = abs(X);P5 = sum(X0.*X0)/N
> t=19.748 sum(X.*conj(X))/N
> t= 3.115 X'*X/N
César Delgado González
Answering first post:
FFT preserve energy, but you should multiply time energy for
(length(Signal)).^2.
For calculating energy signal on time domain (variance), exists
two ways:
Let me call N = length(Signal)
1.- Variance as Sum of samples^2:
N.*sum(Signal.^2)
2.- Variance as std^2:
(N.*std(Signal)).^2
First path is more acurate than the second.
Try this code, see how it works and compare the two ways:
Many speed improvements could be achieve using other calculation
techniques.
Hope this helps,
César