Astrophysicist at work

Just in case I forget how to use google here’s a good page on how to compress eps files – useful in particular when having to submit things on arXiv.

matthew random

I’ve created a Matlab function that properly normalises a histogram with an upper and lower boundary. The main reason for this is that if you have a relatively small number of bins and try to normalise using the trapezium rule just using the area under those bins then you’ll often slightly underestimate the area. This can also be particularly prevalent if your histogram is up against a hard bound, but with a large value e.g. imagine a histogram of positive data peaked close to zero. My function just adds bits on the edges of the histogram data to try and mitigate this effect a bit. Here it is (hopefully self-explanatory):

function [n, x] = histnormbounds(y, m, low, high)

% [n, x] = histnormbounds(y, m, low, high)
%
% Create a normalised historam (normalised using the trapezium rule to
% calculated the area). The low and high values are the boundaries of the
% histogrammed data, and ensure that the histogram gets properly normalised
% if close to the edges of the boundaries. m is the number of histogram
% bins or a vector of bins values.

% set default bounds if none are set (low = -inf and high = inf)
if ~exist('low', 'var')
    low = -inf;
end

if ~exist('high', 'var')
    high = inf;
end

% check that data is within bounds
if min(y) < low && max(y) > high
    error('Data is outside of bounds');
end

% get histrogram
if isscalar(m)
    if m < 2
        error('Need more points in histogram');
    end

    [n, x] = hist(y, m);
else
    if length(m) < 2
        error('Need more points in histogram');
    end

    n = hist(y, m);
    x = m;
end

botbinwidth = x(2)-x(1);
topbinwidth = x(end)-x(end-1); % just in case x isn't uniform

% if histogram points are not close to boundaries (i.e. within a bin of the
% boundaries) then add zeros to histrogram edges
if x(1)-botbinwidth > low
    n = [0, n];
    x = [x(1)-botbinwidth, x];
end

if x(end)+topbinwidth < high
    n = [n, 0];
    x = [x, x(end)+topbinwidth];
end

% if the histogram is closer to the boundary edge than the bin width then
% set a new bin on the boundary with a value linearly extrapolated from the
% gradiant of the adjacent points
if x(1)-botbinwidth < low
    dx = x(1) - low;
    x = [low, x];

    dn = n(2)-n(1);

    nbound = n(1) - (dn/botbinwidth)*dx;

    n = [nbound, n];
end

if x(end)+topbinwidth > high
    dx = high - x(end);
    x = [x, high];

    dn = n(end)-n(end-1);

    nbound = n(end) + (dn/topbinwidth)*dx;

    n = [n, nbound];
end

% normalise n by area
n = n/trapz(x, n);

matthew Matlab

matthew bash, scripts, SFTs

Just a note to say I finally found out (via here) how to get Matlab imagesc plot to have the y-axis with values running from smallest to largest from bottom to top i.e. the correct way! The above link says this info is actually included in the imagesc documentation, but only for a slightly more recent version of Matlab than I have at work. I’ve written a little function to do this for me, which I will share:

function h = imagesc_normal(x, y, N)

% h = imagesc_normal(x, y, N)
%
% This function will return an imagesc-style plot, but with the y-axis
% running from smallest number at the bottom to largest number at the top.
% x is a vector of i values for the x-axis bins, y is a vector of j values
% for the y-axis bins, and N is an i x j array of bin values.
%
% An example of its use would be to plot the output of hist3 e.g.
%   [N, xy] = hist3(randn(1000,2), [20 20]);
%   imagesc_normal(xy{1}, xy{2}, N);

% transpose N
N = N.';

% plot image
h = imagesc(x, y, N);

% set y-axis to correct orientation
set(gca, 'YDir', 'normal');

matthew Matlab

A couple of months back I was trying to figure out how to plot some data on a specific map projection in Matlab. It took me far more googling than should have been necessary, but I eventually go there, so here’s a little reminder to me:

If I have a set of sky positions with their right ascension (between 0 and 2π radians) held in ra and declination (between -π/2 and π/2 radians) held in dec then I could plot those points on a Hammer projection via:

% source RA
ra = 0.6;

% source declination
dec = -0.45;

% create figure
figure;

% create sky map using Hammer projection
axesm('MapProjection', 'hammer', 'AngleUnits', 'radians');

% plot point on sky
plotm(dec, ra, '.');

Or, if I have an image array imarr (for declination running over the first index and right ascension over the second), with points defined at right ascension ras and declination decs, then I could plot the image on the Hammer projection via:

axesm('MapProjection', 'hammer', 'AngleUnits', 'radians');
surfacem(decs, ras, imarr);

More information can be found in Matlab with: doc axesm, doc plotm and doc surfacem.

matthew Matlab

Suppose there was some data of length \(N\) that had been split into \(A\) chunks (each of length \(N/A\)) and each of the chunks had been FFT’d. What I wanted to figure out was how could you recombine these individual FFTs to recreate an FFT of the who N-length data stretch. Now there’s a simple way to do this in Matlab (see here, also linked in struck through stuff below), but I wanted to work out how to do it “properly”, which after a lot of Googling I finally figured out – my method generalises the definition of the decimation-in-frequency radix-4 transform (see e.g. page 32 of this) via the method used for a sequence of two FFTs in this paper.

So, if we start with the definition of the Discrete Fourier Transform and its inverse:
$$X_k = \sum_{n=0}^{N-1}x_nW_N^{nk}$$ and
$$x_k = (1/N)\sum_{n=0}^{N-1}X_nW_N^{-nk}$$
where \(W_N = e^{-2\pi i /N}\).

The decimation-in-frequency FFT for an arbitrary radix (e.g. radix-A) can be defined as:
$$X_{Ak+m} = \sum_{n=0}^{M-1}\left[\sum_{l=0}^{A-1}(W_A)^{lm}x_{n+lM}W_N^{nm}\right]W_M^{nk},$$
where \(M = N/A\), \(k = 0, 1, …, M-1\) and \(m=0,1,…,A-1\). The \(x_{n+lM}\) vectors are the sequential chunks of the original time series. Since we have the FFTs of these chunks we can use the definition of the inverse FFT above and substitute it in for \(x\) giving
$$X_{Ak+m} = \sum_{n=0}^{M-1} \frac{1}{M} \left[\sum_{j=0}^{M-1}\left(\sum_{l=0}^{A-1}W_A^{lm}F_{lj}\right)W_M^{-jn}W_N^{nm}\right]W_M^{nk},$$
where \(F_{lj} = {\rm FFT}(\{x_{jM, jM+1, …, (j+1)M-1}\})_{l}\).

We can simplify this expression with a couple of substitutions, \(W_N = W_M^{1/A}\) and \(d_{jm} = \sum_{l=0}^{A-1} W_A^{lm}F_{lj}\), giving
$$X_{Ak+m} = \frac{1}{M}\sum_{n=0}^{M-1}W_M^n\sum_{j=0}^{M-1}d_{jm}W_M^{(m/A)-j+k},$$
which, using the (I assume standard) identity used in here to get rid of one of the summations), gives
$$X_{Ak+m} = \frac{1}{M}\sum_{j=0}^{M-1}d_{jm}\frac{(1-W_M^{((m/A)-j+k)M})}{(1-W_M^{((m/A)-j+k)})}.$$

And, voila, that’s the general solution to the problem – although note that for the \(Ak+0\) values you should just use
$$X_{Ak+0} = d_{j0},$$
(strangely the general solution doesn’t seem to work for these values, but they may be a problem with what I did).

Below’s a short Matlab code that does this if the \(A\) FFTs have been put into an \(M\times A\) array called X:

% get size of X
sx = size(X);
Nbins = sx(1);
Nffts = sx(2);

% get the length of the final FFT
Ntot = Nbins*Nffts;

% create output Y
Y = zeros(Ntot,1);

m2pi = -2*pi*1i;
tw = exp(m2pi/Nffts);
W = exp(m2pi/Nbins);

jbins = (0:(Nbins-1))';

% get the i*Nbins + 1 full FFT values (for i=0, Nffts)
for i=1:Nffts
    Y(1:Nffts:end) = Y(1:Nffts:end) + X(:,i);
end

% get the rest
for j=1:Nffts-1
    for k=0:Nbins-1
        % get sums of ffts
        ds = zeros(Nbins,1);

        for i=0:Nffts-1
            ds = ds + tw^(j*i)*X(:,i+1);
        end

        % get extra twiddles
        T = (1 - W.^(((j/Nffts) - jbins + k)*Nbins)) ./ ...
            (1 - W.^((j/Nffts) - jbins + k));

        Y(j+1 + Nffts*k) = (1/Nbins)*sum(ds.*T);
    end
end

There are probably quicker ways of implementing this in Matlab. In reality I expect this will always be slower (it seems to be \(\propto N^2/A\) operations) than just inverse FFTing each of the individual FFTs and recombining them to create the full spectrum!

Today I’ve been searching the web trying to find out how you can combine multiple FFT of short chunks of data into an equivalent FFT of the whole stretch of data (e.g. if you have some data that’s 128 point long [let's keep it in powers of two!] and I split it into 4 chunks of 32 points and FFT each of those, how do I combine those FFTs to make the equivalent of an FFT of the whole 128 point long stretch?). After learning a lot about how FFT’s actually work (radixes, butterflies, twiddle factors and all that) I found a nice simple Matlab-fied answer here (See Greg Heath’s reply at – 2010-06-18 22:54:00).

[Update: The method in the link is actually a bit of a cheat as it assumes that the short FFTs have been FFT'd with zero-padding to the full sample length. In reality (in the situation that I want this method for) the short FFTs will not have been zero padded. So, to duplicate the effect of zero padding the short FFTs need to be interpolated with a sinc function to the full sample rate. In Matlab this can be done with the interpft function, although what this really does is (i)FFT the data zero pad it and the re-FFTs it (it's therefore quick), but information on how to do this more formally (using the sinc interpolation) can be found here.]

matthew computing, Matlab

In case I forget here’s how to kill zombie processes (of which many seemed to be spawned, in particular by Chromium, on my Ubuntu 11.04 desktop):

kill -HUP `ps -A -ostat,ppid,pid,cmd | grep -e '^[Zz]' | awk '{print $2}'`

I’ve added this as an alias called killzombie. (Via here, plus multiple other sources).

matthew computing

This is just a reminder to myself that I’ve written up some notes on automating the time domain known pulsar search code on the CW wiki here.

matthew heterodyne code, known pulsar search

After a couple of years of quite disliking the version control software git I’ve finally got used to it and think it’s alright. In fact I’m now using git to access an svn repository with git-svn. Here are just some reminder tips for myself (in part taken from here):

Cloning an svn repository:

git svn clone url_of_svn_repository

Adding and committing to your (local) git repository (this is just like standard git):

git add my_file
git commit -a -m "My commit message"

Committing those commits to the svn repository (like git push):

git svn rebase
git svn dcommit

Removing files from the svn repository (pretty much as above, but the file most be removed first!):

# remove the file
git rm my_file
# commit the removal
git commit -a -m "Removed file"
# commit the changes to the svn repository
git svn rebase
git svn dcommit

There probably a lot more I can do (like adding and merging branches), but this is enough for me at the moment.

matthew computing

As a reminder to myself here’s a bit of LaTeX information on adding comments (that can be turned on/off) in a document. First off use the comment package:

\usepackage{comment}

For comments that enclose stuff that won’t be seen just use the environment:

\begin{comment}
Blah, blah, blah...
\end{comment}

However, if you want to turn comments on or off (so they appear or don’t appear in your final document) you have to place:

\includecomment{comment}, or \excludecomment{comment}

before the \begin{document} command. In fact the “comment” environment name can be anything you want.

You can also have special comments that allow you to change the style of the text within the comment (this is very similar to creating a new environment). So, if I wanted to create a comment environment called notes inside which I wanted all the text to be sans serif and the colour blue I could define:

\specialcomment{notes}{\begingroup\sffamily\color{blue}}{\endgroup}

and within the body of the document use:

\begin{notes}
Here's some stuff that I want to note...
\end{note}

Now, I might want to be able to specify whether to turn the notes on or off when I compile the document. I can do this using the ifthen package and some command line LaTeX e.g. within the document have something like:

\usepackage{comment}
\usepackage{ifthen}
\usepackage{color}

\specialcomment{notes}{\begingroup\sffamily\color{blue}}{\endgroup}

\ifthenelse{\isundefined{\withnotes}}{
  % exclude the comments
  \excludecomment{notes}
}{
  % include the comments
  \includecomment{notes}
}

In this document the \withnotes command is undefined, so by default the notes will be excluded. However, when compiling the document (called, say, mydoc.tex) you can chose to define this command and show the notes e.g.

:~> pdflatex '\newcommand{\withnotes}{}\input{mydoc.tex}'
:~> pdflatex '\newcommand{\withnotes}{}\input{mydoc.tex}'

Or, maybe put this in a Makefile e.g.:

# Make file for mydoc.tex

# default to make document with notes
default: withnotes

withnotes: mydoc.tex
    pdflatex '\newcommand{\withnotes}{}\input{mydoc.tex}'
    pdflatex '\newcommand{\withnotes}{}\input{mydoc.tex}'

withoutnotes: mydoc.tex
    pdflatex mydoc.tex
    pdflatex mydoc.tex

A similar way to do this, but without the easy ability to turn the comments on and off would be to create a new environment like:

\newenvironment{notes}{\sffamily\color{blue}}

matthew LaTeX