Comparison of the derived power spectra and density PDFs for lognormal power distributions. Obviously (if you think about it), the power spectra are always flat for randomly distributed signal. The PDFs change significantly.
In [29]:
img = np.zeros([128,128])
norm = np.random.randn(128,128)
import agpy
lognorm = 10**(norm)*0.1 # peak amplitude ~ 0.1 Jy
lognorm_noisy = lognorm + np.random.randn(128,128)*0.03 # add noise at the 0.03 Jy level (like Bolocam)
clf();
subplot(121); imshow(log10(lognorm)); colorbar()
subplot(122); imshow(log10(lognorm_noisy)); colorbar();
In [30]:
rr,zz = agpy.psds.power_spectrum(lognorm)
rrn,zzn = agpy.psds.power_spectrum(lognorm_noisy)
In [31]:
clf()
subplot(121); loglog(rr,zz,label="Lognormal PSD"); loglog(rrn,zzn,label="Noisy Lognormal PSD")
subplot(122); hist(lognorm.flat,bins=logspace(-3,3)); gca().set_xscale('log')
hist(lognorm_noisy.flat,bins=logspace(-3,3),color='r',alpha=0.5,histtype='stepfilled');
figure();
semilogy(rr,zz,label="Lognormal PSD"); semilogy(rrn,zzn,label="Noisy Lognormal PSD")
Out [31]:
So apparently there's no difference in the power spectrum and only a subtle difference in the PSD. I'll do the same thing with a lower peak, sub-noise
In [33]:
lognorm = 10**(norm)*0.01 # peak amplitude ~ 0.01 Jy
lognorm_noisy = lognorm + np.random.randn(128,128)*0.03 # add noise at the 0.03 Jy level (like Bolocam)
figure();
subplot(121); imshow(log10(lognorm)); colorbar()
subplot(122); imshow(log10(lognorm_noisy)); colorbar();
rr,zz = agpy.psds.power_spectrum(lognorm)
rrn,zzn = agpy.psds.power_spectrum(lognorm_noisy)
figure()
subplot(121); loglog(rr,zz,label="Lognormal PSD"); loglog(rrn,zzn,label="Noisy Lognormal PSD")
subplot(122); hist(lognorm.flat,bins=logspace(-3,3)); gca().set_xscale('log')
hist(lognorm_noisy.flat,bins=logspace(-3,3),color='r',alpha=0.5,histtype='stepfilled');
figure();
semilogy(rr,zz,label="Lognormal PSD"); semilogy(rrn,zzn,label="Noisy Lognormal PSD")
figure()
plot(rr,zz-zzn)
Out [33]:
Any power spectrum differences are exceedingly subtle, at the $10^{-2}$ level
In [34]:
lognorm = 10**(norm)*0.001 # peak amplitude ~ 0.001 Jy
lognorm_noisy = lognorm + np.random.randn(128,128)*0.03 # add noise at the 0.03 Jy level (like Bolocam)
figure();
subplot(121); imshow(log10(lognorm)); colorbar()
subplot(122); imshow(log10(lognorm_noisy)); colorbar();
rr,zz = agpy.psds.power_spectrum(lognorm)
rrn,zzn = agpy.psds.power_spectrum(lognorm_noisy)
figure()
subplot(121); loglog(rr,zz,label="Lognormal PSD"); loglog(rrn,zzn,label="Noisy Lognormal PSD")
subplot(122); hist(lognorm.flat,bins=logspace(-3,3)); gca().set_xscale('log')
hist(lognorm_noisy.flat,bins=logspace(-3,3),color='r',alpha=0.5,histtype='stepfilled');
figure();
semilogy(rr,zz,label="Lognormal PSD"); semilogy(rrn,zzn,label="Noisy Lognormal PSD")
figure()
plot(rr,zz-zzn)
Out [34]:
In [35]:
lognorm = 10**(norm)*1.0 # peak amplitude ~ 1.0 Jy
lognorm_noisy = lognorm + np.random.randn(128,128)*0.03 # add noise at the 0.03 Jy level (like Bolocam)
figure();
subplot(121); imshow(log10(lognorm)); colorbar()
subplot(122); imshow(log10(lognorm_noisy)); colorbar();
rr,zz = agpy.psds.power_spectrum(lognorm)
rrn,zzn = agpy.psds.power_spectrum(lognorm_noisy)
figure()
subplot(121); loglog(rr,zz,label="Lognormal PSD"); loglog(rrn,zzn,label="Noisy Lognormal PSD")
subplot(122); hist(lognorm.flat,bins=logspace(-3,3)); gca().set_xscale('log')
hist(lognorm_noisy.flat,bins=logspace(-3,3),color='r',alpha=0.5,histtype='stepfilled');
figure();
semilogy(rr,zz,label="Lognormal PSD"); semilogy(rrn,zzn,label="Noisy Lognormal PSD")
figure()
plot(rr,zz-zzn)
Out [35]:
No comments:
Post a Comment