Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

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


Math Forum » Discussions » Software » comp.soft-sys.matlab

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

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Freda

Posts: 8
Registered: 1/7/09
Parseval's Theorem
Posted: Jan 10, 2009 11:13 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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.



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.