Hi
I have made a GNU Radio based system to measure the spectrum of the hydrogen
21cm radiation. It is made for RTL-SDR, but also other receivers can be used,
if they have support for GNU Radio.
The main properties of the system are:
* one can plot, in real time, a single spectrum, moving average of a number of
spectra, or an integrated spectrum
* there is also a waterfall plot
* the x-axis of the plot can be frequency, pure Doppler velocity, or Local
Standard of Rest velocity
* one can plot the spectrum into file
* one can write the spectrum into file
* one can make baseband calibration to correct the curvature of the baseband of
a receiver
* one can print the mean value of a spectrum, which can then be used to make
hot/cold calibration to calculate system temperature and scale factor. Scale
factor converts the intensity of the signal from the receiver into antenna
temperature units.
If you want to have accurate intensity units, you want to determine scale
factor which is valid for your setup, because the scale factor depends on the
setup
* one can change the embedded Python block to add other features
Moving average is a nice feature, because the spectrum is updated all the time,
and one can then easily point the antenna to the Galactic plane just by manually
moving an antenna. Of course, the 21cm line is rather noisy, but good enough
to find the Galactic plane.
I like the GNU Radio based approach because one can see everything streamed in
realtime, not just the final, integrated spectrum.
It is fascinating to see the line coming visible from the noise.
I have used this system in GNU Radio version 3.10.7.0, where the RTL_SDR was
controlled via 'osmocom RTL-SDR source', and in GNU Radio version 3.10.12.0,
where the RTL_SDR was controlled via 'Soapy RTLSDR source'. I have noticed that
changing the frequency by using a slider in the flowgraph is much smoother and
faster in the former case.
The Local Standard of Rest velocity calculation gives results which are
consistent with velocity calculators online. However, I was not sure whether the
LSR velocity has to added to or subtracted from the Doppler velocity :-)
Finally, I concluded that it is added. However, if the LSR velocities of the
observed lines are not correct, change the LSR correction from addition to
subtraction.
Attached in this email are the following files:
1) hydrogen_21cm_spectrometer_1024_channels.grc
The flowgraph for GNU Radio. It has 1024 channels in the spectrum.
Note that you cannot change the number of channels in the spectrum just by
changing the parameter vec_length of the flowgraph (see the discussion at
https://github.com/gnuradio/gnuradio/issues/4196)
If you want to have a different number of channels, you have to make
a new version of the flowgraph, located in another directory. See the more
detailed instructions at the end of this email.
I have also implemented the solution given in the GitHub discussion. However,
because it makes the embedded Python block more complicated, I have not included
it here. I will publish it in another email.
2) baseband_spectrum_*.txt
where * is 512, 1024, 2048 or 4096. These files are used when we make the shape
of the baseband flat, because the intrinsic shape of the baseband of RTL-SDR is
quite curved. You need to determine a new baseband spectrum always when you start
a new observing session, because the shape is not constant.
The attached GNU Radio flowgraph needs only the file baseband_spectrum_1024.txt,
the other files are needed only if you change the number of channels.
3) HI-spectrum-plot.pdf
an example of a spectrum plotted by the flowgraph
4) baseband-fit.pdf
an example of a fit of the baseband
5) GNURadio-screenshot-moving-average-spectrum.png
The plot shows moving average over 700 spectra.
The antenna has been manually pointed to the Galactic plane.
6) GNURadio-screenshot-integrated-spectrum.png
An example of an integrated spectrum.
Integration time was about 30 seconds, and the number of spectra integrated was
about 76000.
7) The file functions.py, including the Python function to calculate the Local
Standard of Rest velocity. You also need to have Astropy and Scipy packages
installed in your computer.
All the examples above have been made with a grid-style parabolic antenna with a
size of about 60x100 cm, Nooelec SAWbird+ H1 amplifier, and RTL-SDR V3.
Enjoy the flowgraph :-)
Cheers, Kimmo Lehtinen
-------------------------------------------------------------------------------
Finally, I include here the documentation from the Documentation tab of the
'Process control' block:
1) General
This is a spectrometer for observing the hydrogen 21cm radiation. It has
been designed for the RTL-SDR, but it can be used also with another
receivers, as long as they have GNU Radio support.
This embedded Python block is of type basic block. Thus, there is no
relation between the number of input items (spectra) and the number of
output items. The block has two inputs. One is for single spectra, and one
is for spectra which have gone through the 'moving average' block. The
input items are lists of lists, which means that each item consists of
several spectra, while the number of spectra is varying, typically between
1 and 8 spectra in my system.
The user can choose what to plot in the vector and time raster sinks:
1) single spectra, 2) moving average, or 3) integrated spectrum.
In the cases 1 and 2, only one spectrum of the input items is plotted for
each run of the general_work of this Python block. There is no need to plot
every spectrum, because the stream of data is so large. In the case of 3
above, each spectrum of the input item is used to calculate the averaged
spectrum. Also in this case, only one spectrum is plotted for each run of
the general_work. When one presses the 'Stop integration' button, this
Python block receives all the input items as usual, but they are just
discarded. Instead, the last integrated spectrum is plotted over and over,
until one resets the integration, or goes into other observing mode.
The x-axis of the plot can be 1) pure Doppler velocity, 2) Doppler velocity
corrected for the Local Standard of Rest Velocity (V_LSR) and 3) frequency.
The variables observer_longitude, observer_latitude, galactic_longitude and
galactic_latitude are needed to calculate the right value for the V_LSR.
Of course, the clock of the PC has to be at a right time.
The file functions.py includes the function used to calculate V_LSR. The
LSR velocity is calculated when the flowgraph is started, but it is not
updated automatically after that, because the change in LSR velocity is so
small during a typical integration.
When the flowgraph is running, and you point an antenna to another direction
on the sky, you can change the values of galactic_longitude and
galactic_latitude in the Parameters tab. The LSR velocity is then updated.
However, do not change the values of galactic_longitude and
galactic_latitude when you are integrating.
When starting GNU Radio, the Import block importing the functions from the
file functions.py may give you an error that the Python file cannot be
found. In that case, edit your ~/.bashrc file (or equivalent file in your
system) to include the following two lines:
PYTHONPATH=$PYTHONPATH:/path-to-directory-containg-the-file-functions.py
export PYTHONPATH
Then give a command 'source ~/.bashrc'
The last line in this Python block, 'return 1', tells to GNU Radio that the
output of this Python block consists of only one spectrum. In addition to
spectral data, this Python block also outputs two numbers through the 'QT
GUI Number Sink' blocks.
2) Calculating system noise temperature and conversion factor of the
instrument.
These factors are different for different setups, therefore you want to
derive values which are valid for your setup. To calculate these
parameters, we need to make a so called hot/cold calibration. The
conversion factor scales the raw signal intensity (volts) from the receiver
into antenna temperature units (Kelvin). Follow the steps below: a) Point
the antenna to ground, and integrate for about 10 seconds. Press a button
to print the median value of the spectrum, and note the value which is in
raw units. We call it U_hot in the following. Note the outside temperature
T_amb (ambient temperature, Kelvins). b) Point the antenna to the sky,
away from the Galactic plane, or preferably change the frequency so that
the hydrogen line is outside of the spectrum. Integrate for about 10
seconds, and make sure that neither the hydrogen line, nor any strong
interference signal is visible in the spectrum. Press a button to print the
median value of the spectrum, and note the value which is in raw units. We
call it U_cold in the following. For a brightness temperature of a clear
sky at 1.4 GHz we can use a value of T_sky=15 K.
Calculate the y_factor = U_hot / U_cold
Calculate the system noise temperature
T_sys = (T_amb - y_factor * T_sky) / (y_factor - 1.0)
Calculate the conversion factor G = (U_hot - U_cold) / (T_amb - T_sky)
As an example, I get the values T_sys=379 K and G=1.37e-6.
Update the conversion factor to the variable 'conversion_factor' in the
flowgraph of the GNU Radio. The system noise temperature is not used in
this flowgraph.
3) Correcting the baseband shape of the receiver.
The RTL-SDR receiver has a rather curved baseband. It has to be made flat,
if we want to detect a weak and wide signal, like the hydrogen 21cm
spectrum. The shape of the baseband is given in the accompanied file
baseband_spectrum_1024.txt. The problem with the shape of the baseband is
that it is not constant. For example, the shape changes when one changes
the RF gain of RTL-SDR. Therefore, you want to frequently determine the
shape. Follow the steps below: Point the antenna to the ground, or to the
sky away from the Galactic plane, or change the frequency so that the
hydrogen line is just outside of the spectrum. You may test which of these
methods gives the flattest baseline. Integrate for about 30 seconds, and
make sure that neither the hydrogen line, nor any strong interference
signal is visible in the spectrum. Press the button 'Baseband correction'.
This embedded Python program then derives a calibration spectrum by fitting
the observed spectrum with a Savitzky-Golay function. We want to fit the
observed spectrum with a smooth function, so that there is no noise in the
calibration spectrum. The Python program then scales the fitted function
to a mean value of unity. When we divide the observed spectrum with this
calibration spectrum, the observed spectrum will be flat. However, the very
ends of the spectrum may still be curved. The Python program also plots
the observed spectrum, together with the fit and residuals, to a file with
a name 'baseband_fit.pdf' (or some other file format, as defined by the
user). Note that the shape of the observed spectrum in the plot is much
more curved than usually, because it is without baseband correction, that
is, as it looks coming directly from the RTL-SDR. Check that the fit is
good. The Pythom program removes noise spikes before the fit, but there may
still be noise spikes which affect the fit, particularly if the spikes are
wide. The Python block will apply the new baseband correction
automatically after the fit.
4) About RF gain values
Depends on the specific RF front end chip. The Rafael Micro-based chips
(R820T, R820T2, R828D) have allowable RF gain values of 0.0, 0.9, 1.4, 2.7,
3.7, 7.7, 8.7, 12.5, 14.4, 15.7, 16.6, 19.7, 20.7, 22.9, 25.4, 28.0, 29.7,
32.8, 33.8, 36.4, 37.2, 38.6, 40.2, 42.1, 43.4, 43.9, 44.5, 48.0, and 49.6
dB. However, if provided a gain value within the upper and lower limits
that does not match one of these values, the RTL-SDR will typically select
the next allowable lower value
(from
https://wiki.gnuradio.org/index.php?title=RTL-SDR_FM_Receiver)
5) Number of channels in the spectrum
Do not change the parameter 'vec_length', which is the number of channels
in the spectrum. With the current configuration of this Python block, GNU
Radio cannot handle the change. See the discussion at
https://github.com/gnuradio/gnuradio/issues/4196 If you want to change the number of channels, first copy the files
hydrogen_21cm_spectrometer_1024_channels.grc and
baseband_spectrum_*.txt into a new directory. Then open the copied flowgraph
in GNU Radio, and edit the block 'Process control' (double click on it).
Then change the line
def __init__(self, samp_rate=2.4e6, freq=1420.4, vec_length=1024
into
def __init__(self, samp_rate=2.4e6, freq=1420.4, vec_length=2048
or whatever number of channels you want (one of 512, 1024, 2048 or 4096)
Also, change the variable vec_length in the flowgraph to the same value.
6) Notes
* The conversion factor depends on the gain of the receiver. Thus, if you
change the gain, you need to derive a new conversion factor, if you want to
have accurate intensity values.
* If you want a larger frequency range, change the start and stop values in
the 'QT GUI Range' block, which has an ID equal to 'freq'.
* The 'scientific' way to show spectra is to have Doppler velocity in
x-axis, increasing from left to right. Therefore the frequency in the QT
GUI Vector Sink decreases from left to right, so that the shape of the
spectrum is not inverted.
* The waterfall display in the QT GUI Time Raster Sink block shows the
frequency relative to the center frequency, because in this block the
x-axis start and end values cannot be changed during runtime.
* GNU Radio gives an error 'PLL not locked !' after changing the frequency.
You can ignore it.
* When changing the frequency with a slider, GNU Radio may give errors 'O'
and 's', which mean 'overrun' and 'packet loss between host and device'.
You can ignore them.
* By default, the plots of GNU Radio have 'auto scale' on. You can better
see the effects of moving average and integration by setting auto scale off
(use middle mouse button).
* The purpose of the parameter 'Keep 1 in N': sometimes it is fun to see in
'slow motion' how the noise of the spectrum decreases when integrating the
spectrum. It is also useful if you have an old and slow PC. If the
parameter is for example 10, the flowgraph will process only one spectrum
out of 10.
* In the file functions.py I have ignored all the warnings from Astropy, so
that there are no warnings when Astropy tries to download some parameters
from internet, when the computer is not online.