Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Parseval's Theorem

Subject: Parseval's Theorem

From: Freda

Date: 10 Jan, 2009 16:13:02

Message: 1 of 18

HI,
In the following code, i have implemented a Fourier Transform, and then checked that Parseval's theorem (essentially conservation of energy) holds by summing over all matrix elements of the intensity before and after the FT. For some reason this doesn't seem to hold. Anyone know why?

close all
clear all
%% Create grid of random numbers between 0 and pi.
phase_grid=rand(128)*pi;
subplot(2,4,1),imagesc(phase_grid),title('Phase ranging between 0 and\pi')

%% Create grid of random numbers between 0 and 1.
intensity_grid=rand(128);
%%(to check parseval's thm, sum over intensity)
sum1=sum(sum(intensity_grid,1),2)
sq=sqrt(intensity_grid);
subplot(2,4,2),imagesc(intensity_grid),title('Intensity ranging between 0 and 1')
%%Combine phase and intensity into one complex function.
complex_fn=sq.*exp(i.*phase_grid);
%%Fourier transform the single complex function.
FT=fft2(complex_fn);
%Obtain the phase and intensity of the fourier tranform.
int=abs(FT).^2;
%(returning to parseval's thm, sum over intensity of the fourier transform.
%This should be the same as the sum above.)
sum2=sum(sum(int,1),2)
%The zero pixel needs to be scaled.
int(1,1)=0;
phase=angle(FT);

subplot(2,4,3),imagesc(int),title('Fourier Transform of Intensity')
subplot(2,4,4),imagesc(phase),title('Fourier Transform of Phase')

%%We now filter out high frequencies
min=0;
max=8;
NewFT=ones(128);
for x=1:128
   for y=1:128
        r=sqrt(x.^2 + y.^2);
        if r>max
          NewFT(x,y)=0;
        elseif r<=max
            NewFT(x,y)=FT(x,y);
        elseif r<min
            NewFT(x,y)=0;
        end
   end
end
NewFT(1,1)=0;
a=(abs(NewFT).^2);

subplot(2,4,5),imagesc(a),title('Intensity of FT without higher frequency')
% Now perform and inverse FT on the complex function.
IFT=ifft2(NewFT);
%% Now extract and plot the intensity and phase.
Intensity=abs(IFT).^2;
subplot(2,4,6),imagesc(Intensity),title('Smoothed out intensity')
Phase=angle(IFT);
subplot(2,4,7),imagesc(Phase),title('Smoothed out phase')



Also, on the side, how do i convert all the graphs to grayscale? I tried imshow on a particular graph, and it kind of works, but loses some stuff. There must be a neater way.

Thanks.

Subject: Parseval's Theorem

From: Matt

Date: 10 Jan, 2009 16:44:02

Message: 2 of 18

"Freda " <fredawerdiger@hotmail.com> wrote in message <gkahee$1t3$1@fred.mathworks.com>...
> HI,
> In the following code, i have implemented a Fourier Transform, and then checked that Parseval's theorem (essentially conservation of energy) holds by summing over all matrix elements of the intensity before and after the FT. For some reason this doesn't seem to hold. Anyone know why?

MATLAB's fft functions are not normalized so as to be orthonormal operators. Therefore, Parseval's theorem will not hold until you normalize them. To normalize them divide fft(A) by numel(A), as in the following example


>> a=rand(1,10); norm(a(:)).^2

ans =

    4.8181

>> A=fft(a); norm(A(:)).^2/numel(A) %Check Parseval

ans =

    4.8181

Subject: Parseval's Theorem

From: Matt

Date: 10 Jan, 2009 16:50:03

Message: 3 of 18

"Matt " <mjacobson.removethis@xorantech.com> wrote in message <gkaj8i$qjj$1@fred.mathworks.com>...

> MATLAB's fft functions are not normalized so as to be orthonormal operators. Therefore, Parseval's theorem will not hold until you normalize them. To normalize them divide fft(A) by numel(A), as in the following example

Sorry, I meant normalize as fft(A)/sqrt(numel(A))

Also, to change to grayscale do

>> colormap(gray)

See "help colormap" and "help gray" to fine tune this.

Subject: Parseval's Theorem

From: Freda

Date: 11 Jan, 2009 00:23:01

Message: 4 of 18

"Matt " <mjacobson.removethis@xorantech.com> wrote in message <gkajjr$h1q$1@fred.mathworks.com>...
> "Matt " <mjacobson.removethis@xorantech.com> wrote in message <gkaj8i$qjj$1@fred.mathworks.com>...
>
> > MATLAB's fft functions are not normalized so as to be orthonormal operators. Therefore, Parseval's theorem will not hold until you normalize them. To normalize them divide fft(A) by numel(A), as in the following example
>
> Sorry, I meant normalize as fft(A)/sqrt(numel(A))
>
> Also, to change to grayscale do
>
> >> colormap(gray)
>
> See "help colormap" and "help gray" to fine tune this.
>
>
Matt,
Thankyou so much!

Subject: Parseval's Theorem

From: Greg Heath

Date: 11 Jan, 2009 10:49:17

Message: 5 of 18

On Jan 10, 11:13=A0am, "Freda " <fredawerdi...@hotmail.com> wrote:
> HI,
> In the following code, i have implemented a Fourier Transform, and then c=
hecked that Parseval's theorem (essentially conservation of energy) holds b=
y summing over all matrix elements of the intensity before and after the FT=
. For some reason this doesn't seem to hold. Anyone know why?

FYI:

http://groups.google.com/group/comp.soft-sys.matlab/msg/2222327db2ea7f51

Hope this helps.

Greg

Subject: Parseval's Theorem

From: Matt

Date: 12 Mar, 2009 21:45:19

Message: 6 of 18

Just to make sure the dust settled on this issue, I'm putting some examples of how to orthonormalize MATLAB's FFT and IFFT functions (i.e. to make them consistent with Parseval's theorem)

a=rand(10);

A=fftn(a); A=A/sqrt(numel(A));

anew=ifftn(A); anew=anew*sqrt(numel(anew));

>> norm(a(:))^2, norm(A(:))^2, norm(anew(:))^2

ans =

   36.6115


ans =

   36.6115


ans =

   36.6115

Subject: Parseval's Theorem

From: Yuji Zhang

Date: 13 May, 2013 10:55:07

Message: 7 of 18

"Yuji Zhang" <yuji.zhang@tufts.edu> wrote in message <kmqg9f$ffb$1@newscl01ah.mathworks.com>...
> Hi thanks guys for the discussion~
>
> How about if you want to calculate the energy, which is the area under the curve?
>
> suppose:
> time domain: function a(:), sampling increment: dt
> frequency domain: spectrum A(:), sampling increment: df = 1/(N*dt)
>
> Energy = sum(a.^2)*dt = sum(A.^2)*df = sum(A.^2)/(N*dt)
mistake, should be
Energy = sum(abs(a).^2)*dt = sum(abs(A).^2)*df = sum(abs(A).^2)/(N*dt)


>
> should A = ?
>
> a nice FFT should be A = fftshift(fft(fftshift(a))*dt, right?
>
>
> Thank you~
>
> "Matt J" wrote in message <gpbvpf$ivd$1@fred.mathworks.com>...
> > Just to make sure the dust settled on this issue, I'm putting some examples of how to orthonormalize MATLAB's FFT and IFFT functions (i.e. to make them consistent with Parseval's theorem)
> >
> > a=rand(10);
> >
> > A=fftn(a); A=A/sqrt(numel(A));
> >
> > anew=ifftn(A); anew=anew*sqrt(numel(anew));
> >
> > >> norm(a(:))^2, norm(A(:))^2, norm(anew(:))^2
> >
> > ans =
> >
> > 36.6115
> >
> >
> > ans =
> >
> > 36.6115
> >
> >
> > ans =
> >
> > 36.6115

Subject: Parseval's Theorem

From: Yuji Zhang

Date: 13 May, 2013 18:18:10

Message: 8 of 18

I checked this:

g = time-domain function
G = fftshift(fft(fftshift(g)))*dt;
df = 1/(N*dt);
>> sum(g.*conj(g))*dt
ans =
    0.0707
>> sum(G.*conj(G))*df
ans =
    0.0707

seems it's working. Can anyone familiar with this confirm/comment/correct this? Thanks a lot~

Subject: Parseval's Theorem

From: Yuji Zhang

Date: 13 May, 2013 18:21:16

Message: 9 of 18

Hi everyone I tried this:

g = time-domain function
G = fftshift(fft(fftshift(g)))*dt;
df = 1/(N*dt);
>> sum(g.*conj(g))*dt
ans =
    0.0707
>> sum(G.*conj(G))*df
ans =
    0.0707

I got "G = fftshift(fft(fftshift(g)))*dt;" form a book.

I'm not sure if this method is actually equivalent to the scaling method you mentioned:
a = fft(a)/sqrt(numel(a))

Can somebody familiar with this explain? Thank you!

Subject: Parseval's Theorem

From: Albert

Date: 4 Jun, 2013 10:15:11

Message: 10 of 18

Hello everybody,

while I get the orthonormalization method I don't really understand why we do it... I mean... I try to look at the classic definition of the DFT for example and there we would not divide e^(-j*w*t) by sqrt(numel) right?
I guess my question is why I have never had this problem on paper and I only have it now that I am doing real computations... I mean I have always seen parseval's equation in books with no normalization needed... the only difference I see is that if we consider the continuous frequency domain TF we should integrate over it instead of adding, however I tried to do that with the result of fft with traps() and then multiplying by the step wich would be 2pi/numel... but this of course is not the same as it has been said in this thread as I am dividing by numel and not by sort(numel)....

Any insight will be appreciated.

Subject: Parseval's Theorem

From: Yuji Zhang

Date: 4 Jun, 2013 18:18:08

Message: 11 of 18

Hi Albert~

I'm myself learning too. Just wanna discuss and share my ideas.

Basically, sqrt(numel(f)) roots from the relation between dt and df.
Think of the definition of FT:
F = sum (f * exp(-i*2*pi*f*x) dt) (eq.1)
f = sum (F* exp( i*2*pi*f*x) df) (eq.2)
There is df = 1/(N*dt). (eq.3)
N = numel(f)

dt and df are factored out and we have
F = dt * sum (f * exp(-i*2*pi*f*x)) (eq.4)
f = 1/(N*dt) * sum (F* exp( i*2*pi*f*x)) (eq.5)

The Matlab function fft and ifft does NOT do the "dt" part in eq.3 and 4. - I think this is because you type in fft(x) in matlab, you haven't provided "dt" to matlab.
Matlab does do the 1/N part in eq.4, because matlab automatically know N = numel(x) when you type in fft(x).

So eq. 3 and 4 becomes this:
F = dt * fft (f) (eq.6)
f = 1/dt * ifft (F) (eq.7)

if you wanna use df in eq.6, you do
f = N*df* ifft (F) (eq.8)

This is what i read in the book http://books.google.com/books/about/Numerical_Simulation_of_Optical_Wave_Pro.html?id=RqPZSAAACAAJ. Let me know if you want scan of these couple pages.

I think what people in the thread are doing is, they rescale eq. 6 and 8 as:
F = dt * fft (f) / sqrt(N) (eq.6B)
f = df * ifft (F) * sqrt(N) (eq.8B)
(everything's devided by sqrt(N), to have a symmetry.)

Let me know if this makes sense, everybody~

yuji
 

"Albert" wrote in message <kokenf$4ih$1@newscl01ah.mathworks.com>...
> Hello everybody,
>
> while I get the orthonormalization method I don't really understand why we do it... I mean... I try to look at the classic definition of the DFT for example and there we would not divide e^(-j*w*t) by sqrt(numel) right?
> I guess my question is why I have never had this problem on paper and I only have it now that I am doing real computations... I mean I have always seen parseval's equation in books with no normalization needed... the only difference I see is that if we consider the continuous frequency domain TF we should integrate over it instead of adding, however I tried to do that with the result of fft with traps() and then multiplying by the step wich would be 2pi/numel... but this of course is not the same as it has been said in this thread as I am dividing by numel and not by sort(numel)....
>
> Any insight will be appreciated.

Subject: Parseval's Theorem

From: Albert

Date: 4 Jun, 2013 18:59:09

Message: 12 of 18

Hi,

thanks for the contribution, I don't truly get how you come to the conclusion that

> There is df = 1/(N*dt). (eq.3)

I am sure it has to be pretty simple but I just don't see it right now, is it just becaus when you transform you do the mapping from 0 to N (let's say) to -1/2 to 1/2 in the freq. domain?
If so, tis is also what I expected to fix with the step modification in the integration of the fft...
I myself am following "Signals&Systems" by oppenheim, willsky and nawab, if you have access and you go to page 380 yo find:

sum[n=-inf;+inf](x^2) = (1/2pi)*int[2pi](X^2*dw)

where x is the discrete sequence as a function of "n" and X is the discrete-time fourier transform, that, even though I might be wrong, I view the fft as a sampling of the function.
So by integrating an interpolation of the fft and adjusting the step I would expect to be able to get the integral at the right side of the equal, the left side beeing straight forward.

Also, I don't really get the definition you put of the FFT, the FFT is purely discrete, from discrete time, to discrete frequency, so there is no dt as far as I know, the same goes for the IFFT. Your reasoning might be right, I don't know, but that part does not make sense to me.

I have been thinking that it miht be easier to just look at the base signals of the decomposition as e^(-i*w*x)/sqrt(numel) as that is the whole point of orthonormality as I understand it. Since you have too guaratee that:

sum[-inf;+inf](base_signal1*base_signal2) is 0 if the base_signal1 != base_signal2 and is 1 if base_signal1 = base_signal2.... My only problem/concern is why this factor is not more standard...

what do you think?


"Yuji Zhang" <yuji.zhang@tufts.edu> wrote in message <kolb10$1ek$1@newscl01ah.mathworks.com>...
> Hi Albert~
>
> I'm myself learning too. Just wanna discuss and share my ideas.
>
> Basically, sqrt(numel(f)) roots from the relation between dt and df.
> Think of the definition of FT:
> F = sum (f * exp(-i*2*pi*f*x) dt) (eq.1)
> f = sum (F* exp( i*2*pi*f*x) df) (eq.2)
> There is df = 1/(N*dt). (eq.3)
> N = numel(f)
>
> dt and df are factored out and we have
> F = dt * sum (f * exp(-i*2*pi*f*x)) (eq.4)
> f = 1/(N*dt) * sum (F* exp( i*2*pi*f*x)) (eq.5)
>
> The Matlab function fft and ifft does NOT do the "dt" part in eq.3 and 4. - I think this is because you type in fft(x) in matlab, you haven't provided "dt" to matlab.
> Matlab does do the 1/N part in eq.4, because matlab automatically know N = numel(x) when you type in fft(x).
>
> So eq. 3 and 4 becomes this:
> F = dt * fft (f) (eq.6)
> f = 1/dt * ifft (F) (eq.7)
>
> if you wanna use df in eq.6, you do
> f = N*df* ifft (F) (eq.8)
>
> This is what i read in the book http://books.google.com/books/about/Numerical_Simulation_of_Optical_Wave_Pro.html?id=RqPZSAAACAAJ. Let me know if you want scan of these couple pages.
>
> I think what people in the thread are doing is, they rescale eq. 6 and 8 as:
> F = dt * fft (f) / sqrt(N) (eq.6B)
> f = df * ifft (F) * sqrt(N) (eq.8B)
> (everything's devided by sqrt(N), to have a symmetry.)
>
> Let me know if this makes sense, everybody~
>
> yuji
>
>
> "Albert" wrote in message <kokenf$4ih$1@newscl01ah.mathworks.com>...
> > Hello everybody,
> >
> > while I get the orthonormalization method I don't really understand why we do it... I mean... I try to look at the classic definition of the DFT for example and there we would not divide e^(-j*w*t) by sqrt(numel) right?
> > I guess my question is why I have never had this problem on paper and I only have it now that I am doing real computations... I mean I have always seen parseval's equation in books with no normalization needed... the only difference I see is that if we consider the continuous frequency domain TF we should integrate over it instead of adding, however I tried to do that with the result of fft with traps() and then multiplying by the step wich would be 2pi/numel... but this of course is not the same as it has been said in this thread as I am dividing by numel and not by sort(numel)....
> >
> > Any insight will be appreciated.

Subject: Parseval's Theorem

From: Matt J

Date: 4 Jun, 2013 21:36:11

Message: 13 of 18

"Yuji Zhang" <yuji.zhang@tufts.edu> wrote in message <kmraus$bkc$1@newscl01ah.mathworks.com>...
> Hi everyone I tried this:
>
> g = time-domain function
> G = fftshift(fft(fftshift(g)))*dt;
> df = 1/(N*dt);
> >> sum(g.*conj(g))*dt
> ans =
> 0.0707
> >> sum(G.*conj(G))*df
> ans =
> 0.0707
>
> I got "G = fftshift(fft(fftshift(g)))*dt;" form a book.
>
> I'm not sure if this method is actually equivalent to the scaling method you mentioned:
> a = fft(a)/sqrt(numel(a))
==================

Hi Yuji,

The scaling you've done here is equivalent to the scaling I mentioned in the early days of this thread. However, I find it very interesting/surprising that it turns out to be so. Basically, what you've shown is that, even though the DFT is only an approximation to the continuous domain Fourier Transform (cf. Albert's remarks), the discrete approximations to the continuous domain Parseval integrals are equal both in time and frequency.

To see that it is equivalent, though, it is somewhat easier to write your time domain energy expression as

 E_t = norm(g)^2 * dt

and your frequency domain energy expression as

 E_f = norm(ff(t) *dt )^2* df

Using df*dt=1/N, this can be rewritten

 E_f = norm(ff(t))^2 *df*dt^2
          = norm(ff(t))^2 *dt/N
          = norm(ff(t)/sqrt(N) )^2 *dt

Therefore E_t = E_f if and only if

    norm(g) = norm(G/sqrt(N))

which is the scaling rule that I originally gave.

I did not reach it using these considerations, however. I used matrix algebra ideas instead. MATLAB's FFT operation is equivalent to the matrix vector multiplication

  G=F*g

where the matrix F has entries

 F(k,n)= exp(-j*2*pi *(k-1)*(n-1))

It is easy to show that the norm of the columns of this matrix is sqrt(N)

 norm(F(:,k)) = sqrt(sum( F(:,k).* conj(F(:,k)) ))
                   = norm(ones(N,1) )
                   = sqrt(N)

Therefore, you must normalize F by 1/sqrt(N) in order to get orthonormal columns/rows.

Subject: Parseval's Theorem

From: Matt J

Date: 4 Jun, 2013 21:51:09

Message: 14 of 18

EDIT:


Hi Yuji,

The scaling you've done here is equivalent to the scaling I mentioned in the early days of this thread. However, I find it very interesting/surprising that it turns out to be so. Basically, what you've shown is that, even though the DFT is only an approximation to the continuous domain Fourier Transform (cf. Albert's remarks), the discrete approximations to the continuous domain Parseval integrals are equal both in time and frequency.

To see that it is equivalent, though, it is somewhat easier to write your time domain energy expression as

 E_t = norm(g)^2 * dt

and your frequency domain energy expression as

 E_f = norm(ff(t) *dt )^2* df

Using df*dt=1/N, this can be rewritten

 E_f = norm(ff(t))^2 *df*dt^2
          = norm(ff(t))^2 *dt/N
          = norm(ff(t)/sqrt(N) )^2 *dt

Therefore E_t = E_f if and only if

    norm(g) = norm(ff(g)/sqrt(N))

which is the scaling rule that I originally gave.

I did not reach it using these considerations, however. I used matrix algebra ideas instead. MATLAB's FFT operation is equivalent to the matrix vector multiplication

  fft(g) = F*g

where the matrix F has entries

 F(k,n)= exp(-j*2*pi *(k-1)*(n-1))

It is easy to show that the norm of the columns of this matrix is sqrt(N)

 norm(F(:,k)) = sqrt(sum( F(:,k).* conj(F(:,k)) ))
                   = norm(ones(N,1) )
                   = sqrt(N)

Therefore, you must normalize F by 1/sqrt(N) in order to get orthonormal columns/rows. Equivalently fft(g)/sqrt(N) will preserve the norm of g.

Subject: Parseval's Theorem

From: Albert

Date: 5 Jun, 2013 10:16:09

Message: 15 of 18

Hi Matt,

I still don't get how you can use dt and df in expressions that I think are totally discrete both in time and frequency..
I also don't get this

> E_t = norm(g)^2 * dt

do you mean the order 1 norm or order 2? did you maybe mean norm(g.^2)?
basically the sum of the squares right? not the square of the sum

Also this one

> E_f = norm(ff(t) *dt )^2* df

I just dont get it, what is ff? the fourier transform? but how come it's variable is t? and why multiply by dt?

Could you elaborate a bit more please?


"Matt J" wrote in message <kolngd$9t8$1@newscl01ah.mathworks.com>...
> EDIT:
>
>
> Hi Yuji,
>
> The scaling you've done here is equivalent to the scaling I mentioned in the early days of this thread. However, I find it very interesting/surprising that it turns out to be so. Basically, what you've shown is that, even though the DFT is only an approximation to the continuous domain Fourier Transform (cf. Albert's remarks), the discrete approximations to the continuous domain Parseval integrals are equal both in time and frequency.
>
> To see that it is equivalent, though, it is somewhat easier to write your time domain energy expression as
>
> E_t = norm(g)^2 * dt
>
> and your frequency domain energy expression as
>
> E_f = norm(ff(t) *dt )^2* df
>
> Using df*dt=1/N, this can be rewritten
>
> E_f = norm(ff(t))^2 *df*dt^2
> = norm(ff(t))^2 *dt/N
> = norm(ff(t)/sqrt(N) )^2 *dt
>
> Therefore E_t = E_f if and only if
>
> norm(g) = norm(ff(g)/sqrt(N))
>
> which is the scaling rule that I originally gave.
>
> I did not reach it using these considerations, however. I used matrix algebra ideas instead. MATLAB's FFT operation is equivalent to the matrix vector multiplication
>
> fft(g) = F*g
>
> where the matrix F has entries
>
> F(k,n)= exp(-j*2*pi *(k-1)*(n-1))
>
> It is easy to show that the norm of the columns of this matrix is sqrt(N)
>
> norm(F(:,k)) = sqrt(sum( F(:,k).* conj(F(:,k)) ))
> = norm(ones(N,1) )
> = sqrt(N)
>
> Therefore, you must normalize F by 1/sqrt(N) in order to get orthonormal columns/rows. Equivalently fft(g)/sqrt(N) will preserve the norm of g.

Subject: Parseval's Theorem

From: Matt J

Date: 5 Jun, 2013 14:09:09

Message: 16 of 18

"Albert" wrote in message <kon359$sm1$1@newscl01ah.mathworks.com>...
> Hi Matt,
>
> I still don't get how you can use dt and df in expressions that I think are totally discrete both in time and frequency..
==========

Multiplying discrete sums by sampling interval lengths (dt or df) is a way of numerically approximating a continuous domain integral. For example

 dt=0.1;
 sum(exp(0:dt:1))*dt

is a numerical approximation of the integral of the continuous function exp(t) from t=0 to 1. Similarly,

 sum(abs(g).^2)*dt

is a numerical approximation to the integral of abs(g(t))^2, which is the integral relevant to Parseval's Theorem in continuous time.

> I also don't get this
>
> > E_t = norm(g)^2 * dt
>
> do you mean the order 1 norm or order 2? did you maybe mean norm(g.^2)?
> basically the sum of the squares right? not the square of the sum
===============

MATLAB has a norm() command. I intended wherever I wrote 'norm(g)' to be interpreted as what you would get if you invoked that command on g. By default norm() gives the l2 norm so, yes, it is equivalent to sum(abs(g).^2).



> Also this one
>
> > E_f = norm(ff(t) *dt )^2* df
>
> I just dont get it, what is ff? the fourier transform? but how come it's variable is t? and why multiply by dt?
==============

Sigh. A typo. It should read

   E_f = norm(fft(g) *dt )^2* df

Yuji defined G=fft(g)*dt which, similar to what I mention above, is the discrete numerical approximation to the continuous-time Fourier transform integral. I simply substituted Yuji's original definition in for G. So, going through the derivation again

 E_f= norm(fft(g) *dt )^2* df
     = norm(fft(g))^2*dt^2*df
     = norm(fft(g))^2*dt/N
     = norm(fft(g)/sqrt(N))^2*dt

Equating this to E_t leads to

   norm(fft(g)/sqrt(N)) =norm(g)

Subject: Parseval's Theorem

From: Yuji Zhang

Date: 7 Jun, 2013 00:57:15

Message: 17 of 18

Hi Albert~

I check that book briefly. I don't understand the frequency there is omega = 0 to 2pi.

I can introduce why df = 1/(N*dt)
I think it roots from Nyquist sampling theorem.
http://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem
It says, if your frequency is -B~B Hertz, then your temporal sampling rate should have an increment of 1/2B second.

In a discrete context, we suppose the temporal points are:
linspace(-T/2, T/2, N), the increment is dt = T/N.
Then the frequency we can cover is -1/(2*dt) ~ 1/(2*dt)
( 1/2B = dt => B = 1/(2*dt) )
The frequency length is 1/dt. There are N discrete sample points, so the frequency increment df = 1/(N*dt).



In calculus context, dt, df should be infinitesimal; in discrete context, delta_t, delta_f has certain values. I don't think it matters too much to still use the notations dt, df in the discrete context. You can use delta_t and delta_f to avoid misunderstanding.


Let me think of the other questions. I need to read the book~





"Albert" wrote in message <kolddt$965$1@newscl01ah.mathworks.com>...
> Hi,
>
> thanks for the contribution, I don't truly get how you come to the conclusion that
>
> > There is df = 1/(N*dt). (eq.3)
>
> I am sure it has to be pretty simple but I just don't see it right now, is it just becaus when you transform you do the mapping from 0 to N (let's say) to -1/2 to 1/2 in the freq. domain?
> If so, tis is also what I expected to fix with the step modification in the integration of the fft...
> I myself am following "Signals&Systems" by oppenheim, willsky and nawab, if you have access and you go to page 380 yo find:
>
> sum[n=-inf;+inf](x^2) = (1/2pi)*int[2pi](X^2*dw)
>
> where x is the discrete sequence as a function of "n" and X is the discrete-time fourier transform, that, even though I might be wrong, I view the fft as a sampling of the function.
> So by integrating an interpolation of the fft and adjusting the step I would expect to be able to get the integral at the right side of the equal, the left side beeing straight forward.
>
> Also, I don't really get the definition you put of the FFT, the FFT is purely discrete, from discrete time, to discrete frequency, so there is no dt as far as I know, the same goes for the IFFT. Your reasoning might be right, I don't know, but that part does not make sense to me.
>
> I have been thinking that it miht be easier to just look at the base signals of the decomposition as e^(-i*w*x)/sqrt(numel) as that is the whole point of orthonormality as I understand it. Since you have too guaratee that:
>
> sum[-inf;+inf](base_signal1*base_signal2) is 0 if the base_signal1 != base_signal2 and is 1 if base_signal1 = base_signal2.... My only problem/concern is why this factor is not more standard...
>
> what do you think?
>
>
> "Yuji Zhang" <yuji.zhang@tufts.edu> wrote in message <kolb10$1ek$1@newscl01ah.mathworks.com>...
> > Hi Albert~
> >
> > I'm myself learning too. Just wanna discuss and share my ideas.
> >
> > Basically, sqrt(numel(f)) roots from the relation between dt and df.
> > Think of the definition of FT:
> > F = sum (f * exp(-i*2*pi*f*x) dt) (eq.1)
> > f = sum (F* exp( i*2*pi*f*x) df) (eq.2)
> > There is df = 1/(N*dt). (eq.3)
> > N = numel(f)
> >
> > dt and df are factored out and we have
> > F = dt * sum (f * exp(-i*2*pi*f*x)) (eq.4)
> > f = 1/(N*dt) * sum (F* exp( i*2*pi*f*x)) (eq.5)
> >
> > The Matlab function fft and ifft does NOT do the "dt" part in eq.3 and 4. - I think this is because you type in fft(x) in matlab, you haven't provided "dt" to matlab.
> > Matlab does do the 1/N part in eq.4, because matlab automatically know N = numel(x) when you type in fft(x).
> >
> > So eq. 3 and 4 becomes this:
> > F = dt * fft (f) (eq.6)
> > f = 1/dt * ifft (F) (eq.7)
> >
> > if you wanna use df in eq.6, you do
> > f = N*df* ifft (F) (eq.8)
> >
> > This is what i read in the book http://books.google.com/books/about/Numerical_Simulation_of_Optical_Wave_Pro.html?id=RqPZSAAACAAJ. Let me know if you want scan of these couple pages.
> >
> > I think what people in the thread are doing is, they rescale eq. 6 and 8 as:
> > F = dt * fft (f) / sqrt(N) (eq.6B)
> > f = df * ifft (F) * sqrt(N) (eq.8B)
> > (everything's devided by sqrt(N), to have a symmetry.)
> >
> > Let me know if this makes sense, everybody~
> >
> > yuji
> >
> >
> > "Albert" wrote in message <kokenf$4ih$1@newscl01ah.mathworks.com>...
> > > Hello everybody,
> > >
> > > while I get the orthonormalization method I don't really understand why we do it... I mean... I try to look at the classic definition of the DFT for example and there we would not divide e^(-j*w*t) by sqrt(numel) right?
> > > I guess my question is why I have never had this problem on paper and I only have it now that I am doing real computations... I mean I have always seen parseval's equation in books with no normalization needed... the only difference I see is that if we consider the continuous frequency domain TF we should integrate over it instead of adding, however I tried to do that with the result of fft with traps() and then multiplying by the step wich would be 2pi/numel... but this of course is not the same as it has been said in this thread as I am dividing by numel and not by sort(numel)....
> > >
> > > Any insight will be appreciated.

Subject: Parseval's Theorem

From: Yuji Zhang

Date: 7 Jun, 2013 02:09:09

Message: 18 of 18

Hi Matt~

Nice derivation! I now get it.

The difference is, in your old post you didn't *dt or *df when you calculate the energy. So you do /N (in fft) or *N (in ifft) to scale it.

So let me compile it. It's either the thing in your old post:
> > a=rand(10);
> >
> > A=fftn(a); A=A/sqrt(numel(A));
> >
> > anew=ifftn(A); anew=anew*sqrt(numel(anew));
> >
> > >> norm(a(:))^2, norm(A(:))^2, norm(anew(:))^2
> >
> > ans =
> >
> > 36.6115
> >
> >
> > ans =
> >
> > 36.6115
> >
> >
> > ans =
> >
> > 36.6115

Or this:

>> g = rand(10,1);
>> dt = 0.1;
>> df = 1/(10*dt);
>> G = fft (g) *dt;
>> gnew = ifft(G)*length(G)*df;
>> norm(g)^2*dt, norm(G)^2*df, norm(gnew)^2*dt

ans =
    0.2958

ans =
    0.2958

ans =
    0.2958




Thanks a lot for your help man~

yuji


"Matt J" wrote in message <kolngd$9t8$1@newscl01ah.mathworks.com>...
> EDIT:
>
>
> Hi Yuji,
>
> The scaling you've done here is equivalent to the scaling I mentioned in the early days of this thread. However, I find it very interesting/surprising that it turns out to be so. Basically, what you've shown is that, even though the DFT is only an approximation to the continuous domain Fourier Transform (cf. Albert's remarks), the discrete approximations to the continuous domain Parseval integrals are equal both in time and frequency.
>
> To see that it is equivalent, though, it is somewhat easier to write your time domain energy expression as
>
> E_t = norm(g)^2 * dt
>
> and your frequency domain energy expression as
>
> E_f = norm(ff(t) *dt )^2* df
>
> Using df*dt=1/N, this can be rewritten
>
> E_f = norm(ff(t))^2 *df*dt^2
> = norm(ff(t))^2 *dt/N
> = norm(ff(t)/sqrt(N) )^2 *dt
>
> Therefore E_t = E_f if and only if
>
> norm(g) = norm(ff(g)/sqrt(N))
>
> which is the scaling rule that I originally gave.
>
> I did not reach it using these considerations, however. I used matrix algebra ideas instead. MATLAB's FFT operation is equivalent to the matrix vector multiplication
>
> fft(g) = F*g
>
> where the matrix F has entries
>
> F(k,n)= exp(-j*2*pi *(k-1)*(n-1))
>
> It is easy to show that the norm of the columns of this matrix is sqrt(N)
>
> norm(F(:,k)) = sqrt(sum( F(:,k).* conj(F(:,k)) ))
> = norm(ones(N,1) )
> = sqrt(N)
>
> Therefore, you must normalize F by 1/sqrt(N) in order to get orthonormal columns/rows. Equivalently fft(g)/sqrt(N) will preserve the norm of g.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us