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.

