slightly incorrect calc of bin-to-range conversion

[jkingslake]

The core part of the code that performs the fft on the chirps calculates a conversion between bin number and meters based on the permittivity of ice and the chirp's frequency ramp rate etc. This is called bin2m in the code.

The current version does this in a simplified way which should be improved.

It computes bin2m here:

xapres/xapres/load.py

Line 935 in 00a565d

Profile.bin2m = self.Header["c0"]/(2.*(T1-T0)*pad*math.sqrt(self.Header["ER_ICE"])*self.Header["K"])

i.e

$$ bin2m = \frac{c}{2(T1-T0)p\sqrt{\epsilon}K} $$

This is derived as follows. The frequency bins that come out of the fft are separated in frequency by

$$ f_b = \frac{f_s}{N_{FFT}} $$

where $f_s$ is the sampling frequency (typically 40000 Hz) and $N_{FFT}$ is the number of samples in the fft.

and range and frequency are related by

$$ R = \frac{f_b c}{2\sqrt{\epsilon} K} $$

This is correct so far (I think), but the issue comes when we define $N_{FFT}$.

The code currently uses $N_{FFT} = f_s P (T1-T0)$, where $P$ is the pad factor, $T0$ is the end time of the data used in the chirp (usually 1 second) and $T0$ is the start time of the data used in the chirp (usually 0 seconds).

Equating $R$ with bin2m, we can combine all three equations together to give the expression for bin2m used in the code.

The problem is that this neglects the fact that the data is trimmed by one or two elements here:

xapres/xapres/load.py

Line 928 in 00a565d

winchirp = np.multiply(chirp[0:Nt],np.blackman(Nt))

It is also a bit of a non-elegant way of doing this. It would be simpler not to substitute in an expression for $N_{FFT}$ at all and just carry $\frac{f_s}{N_{FFT}}$ through to the computation of bin2m.

solution

replace

Profile.bin2m = self.Header["c0"]/(2.*(T1-T0)*pad*math.sqrt(self.Header["ER_ICE"])*self.Header["K"])

with

Profile.bin2m = self.Header["c0"]/(2.*math.sqrt(self.Header["ER_ICE"])*self.Header["K"]) / Nfft / self.Header["dt"]

Note that 1/ self.Header["dt"] should equal $f_s$ and that Nfft is the length of the data taking into account the trimming.

[jkingslake]

After looking in more detail it seems like the way the code currently does this is in fact CORRECT in the usual case when the chirp starts out 40001 elements long. In this case it gets trimmed to 40000 elements (before padding), and the expression the code uses is correct in this case because $N_{FFT} = f_s P (T1 − T0)$ is correct.

It is incorrect if the trimming results in a different pre-padding length.

Note on this can be found here https://github.com/ldeo-glaciology/xapres/blob/82e36a71bda23261cc542e1421d319a7acdb45ca/notebooks/test_notes/custom_fft.ipynb

[jkingslake]

after #51 the code does not use this older fft code by default, so this can be closed.