Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Parseval's Theorem
Replies: 18   Last Post: Jun 6, 2013 10:09 PM

 Messages: [ Previous | Next ]
 Freda Posts: 8 Registered: 1/7/09
Parseval's Theorem
Posted: Jan 10, 2009 11:13 AM

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.

Date Subject Author
1/10/09 Freda
1/10/09 Matt
1/10/09 Matt
1/10/09 Freda
1/11/09 Greg Heath
3/12/09 Matt
5/13/13 Yuji
5/13/13 Yuji
5/13/13 Yuji
5/13/13 Yuji
6/4/13 Matt J
6/4/13 Matt J
6/5/13 Albert
6/5/13 Matt J
6/6/13 Yuji
6/4/13 Albert
6/4/13 Yuji
6/4/13 Albert
6/6/13 Yuji