<center><h2>
Class 6
</center></h2>


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

<ul>

<li>Solution to the homework:

<a href="imfilt.m">imfilt.m</a>

<p>We can make a movie as follows.
Note use of logspace function to generate logarithmically
spaced vector of filter cutoffs.

<pre>
im=sum(imread('football.jpg'),3);

r=logspace(0,log10(max(size(im))),50);
for i=1:length(r)
  imagesc(imfilt(im,'low-cut',r(i)));
  title(sprintf('%f',r(i)));
  m(i)=getframe;
end

% play it once
movie(m,0)
</pre>

<li>In this class we're going to make a sample piece of CMB sky.
First get CMB power spectra from output file of standard
program "cmbfast".

<pre>
dat=load('cmbfast_output.txt');
ps.l=dat(:,1);
ps.T=dat(:,2);

% cmbfast is goofy units
ps.Csl=ps.T*(2.73^2)*(1e6^2);
plot(ps.l,ps.Csl);
</pre>

<p>This produces the familiar looking plot.
However the CMB power spectrum is normally given in l(l+1)C_l/2*pi form
- we want straight C_l.

<pre>
ps.Cl=(ps.Csl*2*pi)./(ps.l.*(ps.l+1));
semilogy(ps.l,ps.Cl);
</pre>

<li>We want an array of normally distributed unit variance numbers
which have a purely real Fourier transform.
Try to write function z=randn_ft_of_real(N_pix) which will
do this.
Need to figure out the symmetries required to make ft purely
real...
You know when you have it because the following will return
a purely real result:

<pre>
x=fft2(fftshift(randn_ft_of_real(4)))
</pre>

<p>Here is my version:

<a href="randn_ft_of_real.m">randn_ft_of_real.m</a>

<li>Let's setup a structure containing all the useful
"axis data" for both image and fourier planes:

<pre>
ad.Field_size=10*pi/180;
ad.N_pix=256;
ad.del_u=1/(ad.Field_size);
ad.u_val=ad.del_u*[-ad.N_pix/2:ad.N_pix/2-1];
[x,y]=meshgrid(ad.u_val,ad.u_val);
ad.u_r=sqrt(x.^2+y.^2);
</pre>

<p>Now make some random number and look at them:

<pre>
z=randn_ft_of_real(ad.N_pix);
imagesc(ad.u_val,ad.u_val,real(z)); axis square
</pre>

<li>Now it turns out that there is a simple correspondence between
multipole expansion ell and Fourier u when u is measured
in inverse radians - l=2pi*u

<pre>
ps.u=ps.l/(2*pi);
</pre>

<p>We want to envelope the amplitude of the Fourier modes
so but CMB spectrum is power.
Power is square of amplitude so we do:

<pre>
amp=sqrt(ps.Cl);
</pre>

<li>Now we want to interpolate the 1D power spectrum onto the
2D Fourier grid.
Matlab's interp functions are very powerful and we should take
a minute to look at them.
The interp1 help gives a good example:

<pre>
x = 0:10; y = sin(x); xi = 0:.25:10;
yi = interp1(x,y,xi); plot(x,y,'o',xi,yi,'.')
</pre>

<p>Note that default is linear interpolation.
Other methods are available but need to be used with care
when there is noise in the data.

<pre>
yi = interp1(x,y,xi,'spline');
plot(x,y,'o',xi,yi,'.')
</pre>

<p>Note that interp2 does 2D interpolation.
However what we want to do here is map the 1D cmb power
spectrum onto the 2D Fourier grid.
We already calculated the radius of each point in the Fourier
plane from the origin so we can do an interp1 on that.

<pre>
e=interp1(ps.u,amp,ad.u_r);
</pre>

<p>Because we have a fine pixel resolution in image
space we have a larger span in fourier space - out to ell
numbers beyond the range of our cmbfast run.
interp1 sets out of range interpolations to NaN so
set them to zero and then look at the result:

<pre>
e(~finite(e))=0;

imagesc(ad.u_val,ad.u_val,e); axis square
ca=caxis; caxis(ca/5);
</pre>

<li>This is the "2d amplitude spectrum".
To impose this amplitude spectrum on the random
fourier modes we just multiply the two together:

<pre>
z=z.*e;

imagesc(ad.u_val,ad.u_val,real(z)); axis square
ca=caxis; caxis(ca/5);
</pre>

<li>Taking the Fourier transform of this we get
a piece of sky which has the power spectrum of the CMB
and random Fourier phases:

<pre>
im=ifft2(fftshift(z));
</pre>

<p>Now calculate the image plane axis ticks:

<pre>
ad.del_t=ad.Field_size/ad.N_pix;
ad.t_val_deg=ad.del_t*[-ad.N_pix/2:ad.N_pix/2-1]*180/pi;
</pre>

<p>And finally look at the result:

<pre>
imagesc(ad.t_val_deg,ad.t_val_deg,im); axis square
</pre>

<li>Homework is to:

<ul>
<li>Figure out how to turn the z array back into a power spectrum -
i.e. do the reverse fft and then accumulate the mean square of
the fourier modes in annular bins.
Hint: we can do this by making 2 calls to the histogram function
hfill introduced a couple of classes back, and making use
of the optional weight argument to that function (see "help hfill").

<li>Then run a Monte Carlo to generate 50 realizations of CMB
sky and their corresponding power spectra.
Convert the power spectra back to (l(l+1)Cl/2pi form and
make an errorbar plot showing the mean and sdev of the realizations
in each bin with the input cmbfast model overplotted.

<li>For extra credit check the literature to find a formula for the
expected error bar size for given sky patch area and compare to your
calculation.
Or derive the formula yourself from first principles.

<li>For extra extra credit inject white pixel noise into the image plane
CMB map (i.e. just do im=im+c*randn(size(im)) where c is some
constant chosen such that the noise is significant but does
not completely swamp the CMB signal).
Then make 50 realizations of that and look at how the output
spectrum is biased versus the input.
Why does it look the way it does?
</ul>

</ul>

<hr>