<center><h2>
Class 5
</center></h2>


<center><img SRC="../../../colline.gif" ></center>

<ul>

<li>Solution to the homework:

<pre>
% get the data
s=wmap_read('wmap_ptsrc_catalog_p2_3yr_v2.txt')

% plotting in galactic latitude survey clearly incomplete below b=20
clf
axesm('MapProjection','Hammer','Grid','on');
i=abs(s.b)>20;
plotm(s.b(i),s.l(i),'.')

% area from Pole to 20 deg lat is 2pi(1-cos(70d)) steradians
% so area above 20 deg lat is 4pi(1-cos(20d))
sa=4*pi*(1-cos(70*pi/180));

% fill into log bins - 20 per decade
[x,n]=hfill(log10(s.s(:,1)),30,0,1.5);
% compute bin lower and upper edges
bwl=x(2)-x(1);
bc=10.^x;
bl=10.^(x-bwl/2);
bu=10.^(x+bwl/2);
bw=bu-bl;

% fit to power law form with extra multiplier to account for bin
% widths and survey area - i.e. we are fitting directly to sources per
% Jansky per steradian
[p,pe]=matmin('logl',[2,-3],[],'dnds_fitfunc',n,bc,sa.*bw)

% calculate upper and lower confidence limits for n obs per bin
[le,ue]=poisserrors(n);

% now scale all n per bin according to area and bin width -
% i.e. convert to flux f
f.c=n./(sa.*bw);
f.l=le./(sa.*bw);
f.u=ue./(sa.*bw);

% plot and overplot fit
clf
plot(bc,f.c,'.');
line([bc;bc],[f.u;f.l],'Color','b');
set(gca,'YScale','log'); set(gca,'XScale','log');
hold on
plot(bc,powlaw(p,bc),'r');
hold off

xlabel('Flux S (Jy)');
ylabel('dN/dS (N per Jy per Sr)');

% compare this plot to fig 13 of
% WMAP1 C.L. Bennett, et al., 2003, ApJS, 148, 97
% and compare the numbers in p to table 4
</pre>

<li>Let's look at how much the error bars differ from sqrt(n) ones:

<pre>
plot(x,n,'.')
set(gca,'YScale','log');
line([x;x],[le;ue],'Color','b');
hold on
plot(x+0.01,n,'r.')
hold off
line([x;x]+0.01,[n+sqrt(n);n-sqrt(n)+1e-9],'color','r')
ylim([1e-2,1e2]);
</pre>

<p>We see as expected that the error bars differ for small n.
Note that proper poisson error bars are close to symmetric
about dot on log scale, whereas the sqrt(n) error bars are
strongly asymmetric.

<li>Moving on to Fourier Transforms let's first look at simple 1D fft:

<p>Make a noise sine wave and plot it:

<pre>
delt=0.1;
maxt=20;
t=delt:delt:maxt;
p=[0,2,0.5,0.7]
y=p(1)+p(2)*sin(2*pi*(p(3)*t+p(4)));
y=y+0.5*randn(size(y));
plot(t,y,'.');
grid
</pre>

<p>Now we do the "fast Fourier transform" or fft on it:

<pre>
z=fft(y)
</pre>

<p>Note that the fft of a purely real input gives complex output.
Plot the real part:

<pre>
plot(real(z),'.')
</pre>

<p>Note that the zero freq is at the start, then +ve freq increasing to
right to most +ve half way along, then most -ve freq inc to right for
second half.
So we see two spikes for our 0.5Hz sinusoid at each end of the plot.
This funny order is what fft algorithms naturally put out.
Matlab provides fftshift function to put the zero freq in the middle
where it belongs:

<pre>
z=fftshift(z);
plot(real(z),'.')
</pre>

<p>We will work always with the zero freq at the center
and treat the need to fftshift as an annoyance which we
always do after/before doing fft operations.

<li>The mantra when working with fft's is "span goes to resolution
and resolution goes to span". i.e. for a wide span in time we have
fine resolution in frequency, and for fast sample rate we measure
a wide span in frequency (up to very high freq).
This makes sense - if you want to capture very fast oscillations
clearly you have to sample very fast.

<p>Let's work out the frequency axis corresponding to out time
axis t:

<pre>
% delta freq is one over span in time
delf=1/(maxt);

% one freq in fft output for each sample in time input
n=length(t);
f=delf*[-n/2:n/2-1];

% plot freq spectrum with freq axis
plot(f,real(z),'.')
</pre>

<p>Note that the zero freq lives in the n/2+1 element.
We see spikes at +/-0.5 Hz as expected.
But what is negative freq anyway?
We already noted pure real input gives complex output,
but if the input is purely real then the complex output
is "conjugate symmetric" about zero freq.
More specifically the first and n/2+1 elements are of
the fft are purely real and the 2:n/2 and n/2+2:n complex
element are (conjugate) mirror images of one another:

<pre>
[z(1) z(n/2+1)]
[z(2:n/2) z(end)]
[z(n/2) z(n/2+2)]
</pre>

<li>Instead of plotting the real and imag parts of the fft output
separately we are often interested in "how much power is
in each Fourier mode".
Power is the square of voltage and in this case the relevant
quantity is z.*conj(z).

<pre>
p=z.*conj(z);
plot(f,p,'.')
% note that real(z).^2+imag(z).^2 = abs(z).^2 = z.*conj(z)
</pre>

<p>The above plot is the "two sided power spectral density" or PSD.
Since the +ve and -ve freq are identical of course we only need
to plot one side:

<pre>
plot(f(n/2+1:end),p(n/2+1:end),'.');
</pre>

<li>Now moving on to 2D FFT's first let's get an image from a file:

<pre>
% read image and sum RGB to form a grayscale image
im=sum(imread('football.jpg'),3);
% force the image square
n=min(size(im));
im=im(1:n,1:n);
% plot it
colormap gray
imagesc(im); axis square
</pre>

<p>Now fft it - note we use fft2 - if we do fft on 2d array we will
get the 1d fft of each column which is not what we want.

<pre>
z=fft2(im);
</pre>

<p>Now we look at the transform and don't see much because the colorscale
has been forced by the large values in the few lowest freq modes.
So we truncated the color scale so we can see the higher freq modes:

<pre>
% look at the transform
imagesc(real(z)); axis square; ca=caxis; caxis(ca/100);

% shift it
z=fftshift(z);

% look again
% this is what 2d fft's typically look like...
imagesc(real(z)); axis square; ca=caxis; caxis(ca/100);
</pre>

<p>No let's low pass filter the image.
The zero order way to do this is to zero the high freq
modes.
First we gen the freq axis (in arbitrary units)

<pre>
f=-n/2:n/2-1;
</pre>

<p>Now we want to construct an array containing the radius
from the center of the fourier plane (i.e. the absolute
value of the frequency for each 2d mode).

<pre>
% make grids of x/y freqs
[fx,fy]=meshgrid(f,f);

% look at them
imagesc(fx);
imagesc(fy);

% make a grid of radius from center
r=sqrt(fx.^2+fy.^2);

% look at it
imagesc(r)
</pre>

<p>Note the very useful Matlab function meshgrid which constructs
2d arrays from 1d vectors. (It's a special canned double version of
repmat.)
We can now apply hard cutoff low pass filter to fourier plane z:

<pre>
z(r>20)=0;

% look at the result in Fourier space - see "coin" in center
imagesc(real(z)); axis square; ca=caxis; caxis(ca/100);

% undo the fft shift
z=fftshift(z);

% look again - note quarter coins in corners
imagesc(real(z)); axis square; ca=caxis; caxis(ca/100);
</pre>

<p>Now we do the inverse fourier transform using ifft and
view the blurred result:

<pre>
imp=ifft2(z);

% view the blurred result in image space
imagesc(imp); axis equal
</pre>

<li>Now let's take a step back and look at what individual
Fourier modes look like in image space.
First we setup and empty grid, then inject a value into
single mode and then view the result:

<pre>
% size of grid
n=256;
% center of fourier plane - zero freq
c=n/2+1;
% fill grid with zeros
z=zeros(n);
% inject a value into a single F mode
z(c+10,c)=1;
% do the inverse transform
im=fft2(fftshift(z));
% view the result
imagesc(im);
</pre>

<p>imagesc barfs because input is complex - but why?
Because we set only the z(c+10,c) mode and left it's symmetric
partner at zero.
We need to do:

<pre>
% set the partner
z(c-10,c)=1;
% do the inverse transform
im=fft2(fftshift(z));
% view the result
imagesc(im);
grid
</pre>

<p>Now we see a "corrugation" - a pattern which varies sinusoidally
in the one direction and is constant in the orthogonal direction.
Because we set a mode on the y-axis in Fourier plane the
pattern varies in the y direction.
Note that this mode varies as the cosine - it is maximum at the center.
Now let's set the mode to complex value 0+i instead:

<pre>
% zero the array
z=zeros(n);

% set the mode to i
z(c+10,c)=complex(0,1);

% ...and its symmetric partner to -i
z(c-10,c)=conj(complex(0,1));

% do the inverse transform
im=fft2(fftshift(z));
% view the result
imagesc(im);
</pre>

<p>Note that this mode is sine wave - it is zero at center.
Each Fourier mode has real part corresponding to cosine
corrugation of given freq and orientation and imaginary
part corresponding to sine corrugation of same freq and
direction.

<li>Homework is to implement function im=imfilt(im,type,r).
This function accepts a 2d image and returns a filtered version of
that image.
type argument is a string which can be "low-cut" to do
a low-pass hard cutoff filter at freq r,
"high-cut" to do a high-pass hard cutoff filter,
or "low-gauss" to apply a Gaussian low pass rolloff to the Fourier
plane where r is the "sigma" of the Gaussian.
r is in arbitrary units of frequency (see above).

<p>Extra credit for making it accept non-square input image.

<p>Extra, extra credit for making an animation using "getframe"
command  in a loop over imfilt with type "low-cut" and increasing r.
Play the the movie using "movie" command.
If you do this last part do it in a separate file which calls
imfilt function rather than adding additional code within imfilt.

</ul>

<hr>