% load cmbfast spectrum from file in useful form
ps=load_cmbfast()
% generate ad structure defining image/fourier grids
ad=calc_ad(10,256);
% interpolate amp to f grid
e=interp1(ps.u,sqrt(ps.Cl),ad.u_r);
e(~finite(e))=0;
% make 50 realizations of the cmb sky and take their power spectra
for i=1:50
% get a random F plane with real ft
z=randn_ft_of_real(ad.N_pix);
% envelope it
z=z.*e;
% go to image plane
im=ifft2(fftshift(z));
% add image plane noise
if(0)
im=im+std(im(:))*randn(size(im));
end
if(i==1)
figure(1)
imagesc(ad.t_val_deg,ad.t_val_deg,im); axis image
end
% go back to f plane
z=fftshift(fft2(im));
% square to get power estimates
z=z.*conj(z);
% accumulate the total and number of samples in annular bins
max_u=2500/(2*pi);
[psr.bcu,tot] =hfill(ad.u_r,50,0,max_u,z);
[psr.bcu,nhits]=hfill(ad.u_r,50,0,max_u);
% take the ratio of these to get the mean
psr.Cl(:,i)=tot./nhits;
end
% convert to conventional CMB form
psr.l=psr.bcu'*2*pi;
l=repmat(psr.l,[1,size(psr.Cl,2)]);
psr.Csl=psr.Cl.*l.*(l+1)/(2*pi);
% plot with the input spectra overplotted
figure(2)
plot(psr.l,psr.Csl,'.',ps.l,ps.Csl,'r')
% Take Csl error from sim
psr.Csles=std(psr.Csl,[],2);
% and mean value
psr.Csl=mean(psr.Csl,2);
% plot as error bars
figure(3)
errorbar(psr.l,psr.Csl,psr.Csles,'.')
hold on
plot(ps.l,ps.Csl,'r')
hold off
dof=2*20+1; x=1:100; plot(chi2pdf(x,dof))
The mean of a chi2 distribution is the dof and the standard deviation is sqrt(2*dof):
xs=chi2rnd(dof,100000,1); mean(xs) std(xs)
So to get the expected error distribution for full sky coverage if the expectation value is 1000 we just normalize:
Cl=1000; plot(x*(Cl/dof),(Cl/dof)*chi2pdf(x,dof))
The std error bar also just gets multiplied by Cl/dof giving us equation 4 in the link above. If we have less than full sky coverage we simply have correspondingly less degrees of freedom and a bigger error bar as equation 8 of the link. The fraction of the sky in our 10x10deg patch is:
fsky=(10*pi/180)^2/(4*pi);
On the other hand we have binned into bins many multipoles wide giving us a boost to the number of degrees of freedom. We can approximate this by multiplying by the bin width. (Strictly we should integrate over the l bin.) So we have:
dl=psr.l(2)-psr.l(1); psr.Cslet=psr.Csl.*sqrt(2./((2*psr.l+1)*fsky*dl)); hold on errorbar(psr.l+10,psr.Csl,psr.Cslet,'m.'); hold off
But why calculate the number of degrees of freedom so elaborately - we already have it - it is just the number of fourier modes in each of our annular bins:
plot(psr.l,nhits,'b',psr.l,2*psr.l*fsky*dl,'r')
% generate ad structure defining image/fourier grids ad=calc_ad(10,256); % get a random F plane with real ft z=randn_ft_of_real(ad.N_pix); % look at it imagesc(real(z)); axis image % go to image plane im=ifft2(fftshift(z)); % construct a mask with soft edge n=50; m=50; x=0.5+0.5*sin(linspace(-pi/2,pi/2,n)); y=ones(1,ad.N_pix); x=[zeros(1,m),x]; y(1:n+m)=x; y(ad.N_pix-n-m+1:end)=fliplr(x); % make 2d m=y'*y; % mask the image im=im.*m; % go back to f plane zp=fftshift(fft2(im)); % compare the two subplot(1,3,1) imagesc(real(z)); axis image subplot(1,3,2) imagesc(m); axis image subplot(1,3,3) imagesc(real(zp)); axis image
We load in some data from a row of 8 DASI pointings and look at what we've got:
load dvia dvia
The u,v elements are the points in the fourier plane where DASI took a sample.
u=dvi(6).u; v=dvi(6).v; vis=dvi(6).vis; plot(u,v,'+'); axis square
At each of these points we have a complex visibility. From these we can construct an image. We want to bin the obs into a fourier grid. First we append the complex conjugate obs:
u=[u;-u]; v=[v,-v]; vis=[vis;conj(vis)];
Make an "axis data" structure for a 5x5deg field with 256 pixels.
ad=calc_ad(5,256);
Now we use the hfill2 function which is 2D analog of the hfill function to bin into a grid. We need to be a little bit careful to make sure the x_tic/y_tic grid output by hfill2 matches the ad.u_val grid exactly.
n=ad.N_pix; l=ad.u_val(1)-ad.del_u/2; h=ad.u_val(end)+ad.del_u/2; % accumulate mean real and imag in each cell [x_tic,y_tic,nhits]=hfill2(u,v,n,l,h,n,l,h); [x_tic,y_tic,r]=hfill2(u,v,n,l,h,n,l,h,real(vis)); [x_tic,y_tic,i]=hfill2(u,v,n,l,h,n,l,h,imag(vis)); z=complex(r,i)./nhits; % force cells where we have no obs to zero z(isnan(z))=0;
Now we ft to the image plane. Note we need to fftshift the image this time since it is a real image of the sky it matters "where the quadrants are". Also we have to force the output real as numerical precision limits mean that it otherwise has a tiny imag component.
im=real(fftshift(ifft2(fftshift(z)))); % plot the image imagesc(ad.t_val_deg,ad.t_val_deg,im); axis xy % astronomy convention has RA inc to the left set(gca,'XDir','reverse')