Calculating flux density of radio telescope

[GhostLoveScore]

Hello

I've made a small program that does FFT, Fast Fourier Transform on a signal from radio receiver (RTL SDR), and parabolic antenna. Of course, the output is amplitude depending on frequency.
I want to use it to detect and measure 1420MHz radiation from space. But I'm not sure on what's the best way to do it. I know the intensity is proportional to amplitude squared, and intensity I= P/A., where I is intensity, P is power received, and A is area of my antenna.
A common unit in radio astronomy seems to be Jansky (flux density) which is proportional to P/(A*bw) where P is power received, A is area of the receiving antenna and bw is bandwidth of your receiving system.

So, I'm a bit confused... My assumption is that I should measure what my antenna is outputting in Janskys. (is that correct?)
That means:

flux density= I/bw=V^2/bw, where V is output from my receiver. But since I'm doing a FFT over a range of frequencies (about 2MHz), what should I do here? How can I calculate flux density from my FFT result? Integrate the amplitudes over the frequency range?? Something else?

 

[tech99]

Summary:: I've done FFT on a signal from my radio receiver, now I want to use the data to calculate radiation intensity from the sky.

Hello

I've made a small program that does FFT, Fast Fourier Transform on a signal from radio receiver (RTL SDR), and parabolic antenna. Of course, the output is amplitude depending on frequency.
I want to use it to detect and measure 1420MHz radiation from space. But I'm not sure on what's the best way to do it. I know the intensity is proportional to amplitude squared, and intensity I= P/A., where I is intensity, P is power received, and A is area of my antenna.
A common unit in radio astronomy seems to be Jansky (flux density) which is proportional to P/(A*bw) where P is power received, A is area of the receiving antenna and bw is bandwidth of your receiving system.

So, I'm a bit confused... My assumption is that I should measure what my antenna is outputting in Janskys. (is that correct?)
That means:

flux density= I/bw=V^2/bw, where V is output from my receiver. But since I'm doing a FFT over a range of frequencies (about 2MHz), what should I do here? How can I calculate flux density from my FFT result? Integrate the amplitudes over the frequency range?? Something else?

I am not sure if the voltage readings at various frequencies which your system produces can be converted to noise power. Imagine that in the analogue world we used a tunable filter to measure each Fourier component. Each of these components would itself be noise and so its power cannot be measured using a voltmeter or by taking peak a amplitude squared.

 

[Ibix]

Why not Fourier transform, filter out everything except (1420$\pm$whatever)MHz band, inverse Fourier transform and then you can calculate the power in your surviving frequencies easily enough. Or am I missing what you are trying to do?

 

[GhostLoveScore]

Why not Fourier transform, filter out everything except (1420$\pm$whatever)MHz band, inverse Fourier transform and then you can calculate the power in your surviving frequencies easily enough. Or am I missing what you are trying to do?

Yes, that's exactly what I'm trying to do.
Then after I do inverse Fourier transform, I should calculate the power by integrating squares of the amplitude over some time interval?

 

[tech99]

I can see what you are trying to do. Add the squares of all the component voltages and then take the square root, so we end up with the RMS value. My problem is that noise has an infinite number of Fourier components. Expressed another way, each component will be noise rather than a sine wave, so how do we find its power?

 

[alan123hk]

Please note that if you want to measure the power of a random signal, it must be much stronger than the background noise, otherwise it cannot be measured. I believe, if I am not mistaken, as long as you can afford the price, a spectrum analyzer with channel power measurement function can accomplish this task. If the price is too expensive, there seems to be a rental service.

Of course, this is just my personal thought, please consult relevant professionals or sales engineers to understand whether it can really meet your measurement requirements.

How to Make Channel Power Measurements : ​

 

[GhostLoveScore]

I can see what you are trying to do. Add the squares of all the component voltages and then take the square root, so we end up with the RMS value. My problem is that noise has an infinite number of Fourier components. Expressed another way, each component will be noise rather than a sine wave, so how do we find its power?

Like Ibix said, could we filter out everything except some frequency interval around the center frequency, then do inverse Fourier transform, and then do the sum of squares of amplitudes? It still wouldn't be single frequency sine wave, but maybe it would be a bit less noisy when we filter out most of 2MHz bandwidth.

How to Make Channel Power Measurements :

"To find the power of this modulated signal we need to find the area under the curve, in other words, we need to integrate." - That was one of my assumptions in my first post, and seems like an easy way to measure the power.

 

[hutchphd]

Why not Fourier transform, filter out everything except (1420$\pm$whatever)MHz band, inverse Fourier transform and then you can calculate the power in your surviving frequencies easily enough. Or am I missing what you are trying to do?

Given Parseval's Theorem, why would you bother to invert the Fourier Transform? Just choose the appropriate band (or bands with normalization) from the FFT.
Is there a method by which you can calibrate this? There are a lot of steps here.

 

[Ibix]

Given Parseval's Theorem, why would you bother to invert the Fourier Transform?

I presumed the OP would want an answer in the time domain.

Perhaps we should ask the OP what the aim of this game is. And what, exactly, has been Fourier transformed.

 

[alan123hk]

"To find the power of this modulated signal we need to find the area under the curve, in other words, we need to integrate." - That was one of my assumptions in my first post, and seems like an easy way to measure the power.

The spectrum analyzer can display the amplitude spectrum and the power spectrum. With the power spectrum, of course, you can integrate any frequency band to find the total power of this frequency band.

For the data converted from the time domain to the frequency domain with FFT, based on the Parseval's theorem, the power of any frequency band can be directly calculated.

https://en.wikipedia.org/wiki/Parseval's_theorem#:~:text=In mathematics, Parseval's theorem usually,the square of its transform.

 

[GhostLoveScore]

I presumed the OP would want an answer in the time domain.

Perhaps we should ask the OP what the aim of this game is. And what, exactly, has been Fourier transformed.

As I said, I'm recording data from RTL SDR, which records in-phase and quadrature data (IQ). I record N samples, then do FFT on them and get frequency domain data. The goal is to somehow calculate the power received from certain bandwidth around 1420MHz. What I'm hoping is to be able to detect when Milky way is passing in front of my antenna and hopefully detect increase in power received.

The spectrum analyzer can display the amplitude spectrum and the power spectrum. With the power spectrum, of course, you can integrate any frequency band to find the total power of this frequency band.

For the data converted from the time domain to the frequency domain with FFT, based on the Parseval's theorem, the power of any frequency band can be directly calculated.

https://en.wikipedia.org/wiki/Parseval's_theorem#:~:text=In mathematics, Parseval's theorem usually,the square of its transform.

So basically that means I can simply integrate the squares in frequency domain and get received power?

 

[tech99]

The video suggests integrating the frequency response by counting squares. As the display is in decibels this will not work, because decibels cannot be added.
Further, although the spectrum analyser says it displays power, how is this likely? The A to D converter is a voltage measuring device. To measure power in any slice of spectrum we need to use a square law detector; we cannot do it using a linear detector (such as an A to D) followed by a correction to obtain RMS, because we do not know the shape of the wave.

 

[Ibix]

So you are hoping to see the received power from hydrogen alpha increase as the Milky Way passes across your field of view. Ok. Do I understand that you've simply recorded a time series of data and Fourier transformed the entire thing? If so, I don’t think that's the right approach. If you filter it as I suggested you will lose the modulation, which is what you are trying to detect.

What a spectrum analyser does is (simplifying somewhat) Fourier transform the last hundred (or whatever) data points, so you get a time series of the average power received in each frequency range over the preceding 100 sample periods. The more samples you include in the window the finer the frequency ranges you can detect, but the more time you have averaged over and the worse your time resolution is.

I think you want to look up "windowed Fourier transform", and listen to people here who appear to know more about it than I do. 😁

 

[GhostLoveScore]

The video suggests integrating the frequency response by counting squares. As the display is in decibels this will not work, because decibels cannot be added.

Would it work if I would choose linear y-axis instead of logarythmic?

Left one is log, right one linear. Just a notice, these graphs are taken at slightly different times, that's why they look a bit different.

1) Do I understand that you've simply recorded a time series of data and Fourier transformed the entire thing?


2)I think you want to look up "windowed Fourier transform", and listen to people here who appear to know more about it than I do. 😁

1) No, I'm doing a transform on 2048 samples, I forgot to write that. Entire thing is 262144 samples if I remember correctly.

2) I'm using Hann window. I've noticed slightly different looks to the graph with and without the window function so I'm guessing it works. How I'm implementing it is, before FFT I'm multiplying sample value with window function (which is dependent on "n", sample number). I'm multiplying both I and Q components, I'm not sure if that's the way to do it, I'm experimenting until I find out how it works :)

And one more thing. These graphs are made by averaging FFT output from 50 different "passes" of 2048 samples each. When not averaged they are very noisy.
And although the graph says "Hydrogen line" this is, of course not it. The graphs are taken at 93MHz center frequency.

 

[alan123hk]

So basically that means I can simply integrate the squares in frequency domain and get received power?

The mathematical definitions of Discrete Fourier Transform and Inverse Discrete Fourier Transform are as follows.

$$ \text {DFT} ~~~~~~~~~~~~X(k) = \sum_{n=0}^{N-1} x(n) e^{\left( \frac {-j2 \pi kn}{N}\right) } ~~~~~~~~~~~~~~~ \text {For} ~~~ 0 \leq k \leq N-1$$
$$ \text {IDFT} ~~~~~~~~~~x(n) = \frac {1} {N}\sum_{n=0}^{N-1} X(k) e^{ \left( \frac {j2 \pi kn}{N} \right) } ~~~~~~~~~~~~ \text {For}~~~ 0 \leq n \leq N-1$$
$$ P=~\frac {1} {N} \sum_{n=0}^{N-1} |x(n)|^2 ~=~ \frac {1} {N^2} \sum_{n=0}^{N-1} |X(k)|^2$$

 

[hutchphd]

Its all right here also:

https://en.wikipedia.org/wiki/Parse...s theorem usually,the square of its transform.

 

[GhostLoveScore]

Thanks everyone.

 

[Baluncore]

These graphs are made by averaging FFT output from 50 different "passes" of 2048 samples each. When not averaged they are very noisy.

With the Hann window, you should overlap your time windows by 50% so that the input data over time is given equal weight.

Power is proportional to the square of the amplitude of each FFT channel. The plot routine may display in dB, but the raw output of an FFT is pairs of orthogonal amplitude components, sine and cosine, from which you calculate channel power = x² + y². Obviously you do not need to compute the square root to get amplitude, before squaring to get the channel power.

Power spectrum accumulation is the only rational way to average your forward transformed data. That is simply the average of all the power spectrums.

 

[GhostLoveScore]

With the Hann window, you should overlap your time windows by 50% so that the input data over time is given equal weight.

Thanks, I didn't know I needed that, and after reading about it, it's really obvious why it's used.

EDIT:
One more maybe stupid question: Do I need to subtract noise level when calculating power? I assume I do, but I'm asking just in case.

 

[Baluncore]

The advantage of power spectrum accumulation is that the signal will rise out of the noise. There is no advantage in subtracting the noise floor, unless you have a dynamic range limitation in your number representation. You might do better with the display by normalising the peak power accumulated and then plotting the noise as it falls with the √(n samples).
Some FFT output channels will randomly give very low noise powers, so immediate subtraction of the noise floor will give them negative power, which has number system implications to the averaging.

The RTL SDR is a hot receiver, so the noise “floor” will be quite high. You may be able to detect the Sun, when there is a solar storm, but there are a limited number of astronomical sources that will climb out of the hot receiver's front end noise. Identify known sources and set your initial antenna declination to sweep past those sources during the sidereal day. A low noise preamp at the dish focus could make a big difference to the number of detectable sources.

 

[GhostLoveScore]

A low noise preamp at the dish focus could make a big difference to the number of detectable sources.

Yes, I plan on using a preamp for hydrogen line.

 

[Baluncore]

What is your dish diameter and your latitude? That has implications to the best antenna declination for initial experiments. I would aim to get a couple of strong sources passing through the fixed beam each day.

You will succeed because you know enough to ask the right questions and avoid the worst pitfalls. If you don't have a copy yet, get a copy of Radio Astronomy by John Kraus.

 

[GhostLoveScore]

What is your dish diameter and your latitude? That has implications to the best antenna declination for initial experiments. I would aim to get a couple of strong sources passing through the fixed beam each day.

You will succeed because you know enough to ask the right questions and avoid the worst pitfalls. If you don't have a copy yet, get a copy of Radio Astronomy by John Kraus.

I don't have a dish yet. I'm reading into building my own, it should be around 2.5m diameter. But don't take my word for it. I haven't read enough about parabolic antennas to give you 100% definitive answer.

45N latitude. No, I don't have that book, I'll go see on Amazon.

 

[Baluncore]

λ = 300 / 1421; ∴ λ = 210 mm.
Given a 2.4 m dish diameter; 2400 / 210 = 11.43 λ diameter.
57° / 11.43 = 5.0° beamwidth.

You can use 1/2" chicken wire for the dish surface. (Holes less than 20 mm).
Dish surface accuracy only needs to be ±½” or ±12 mm.

Latitude 45°N suggests;
Taurus A at dec = 50°
Cygnus A at dec= 40.5°
Cassiopeia A at dec = 58.5°

 

[GhostLoveScore]

λ = 300 / 1421; ∴ λ = 210 mm.
Given a 2.4 m dish diameter; 2400 / 210 = 11.43 λ diameter.
57° / 11.43 = 5.0° beamwidth.

You can use 1/2" chicken wire for the dish surface. (Holes less than 20 mm).
Dish surface accuracy only needs to be ±½” or ±12 mm.

Latitude 45°N suggests;
Taurus A at dec = 50°
Cygnus A at dec= 40.5°
Cassiopeia A at dec = 58.5°

I think I'll check Cassiopeia A first - brightest radio source in the sky (above 1GHz), makes sense to start with it.

 

[Baluncore]

There is good info here on the high flux sources.
https://www.radio2space.com/cassiopeia-a-recorded-with-spider-radio-telescope/

 

[tech99]

When integrating, the Y axis must be in power and not decibels of power.

 

[GhostLoveScore]

When integrating, the Y axis must be in power and not decibels of power.

I know. If you look at my post, post number 14 https://www.physicsforums.com/threa...nsity-of-radio-telescope.1002515/post-6486692, you will see that I'm using logarithmic scale on purpose, so that it's easier to compare with spectrum I'm seeing in SDR programs which use logarithmic scale.

Right graph is what I'm getting when plotting output from FFT in linear y scale.

If anybody is interested, I'm using fftw library to do FFT.

 

[Baluncore]

When integrating, the Y axis must be in power and not decibels of power.

Is power spectrum accumulation really integration?
I say no, the power estimates are discrete in time.

What would happen if you averaged the dB power? If you accumulate dBs you are summing log(power); which is equivalent to multiplication in the linear domain. When you divide that sum by n, you are effectively extracting the n'th root.

The Arithmetic Mean is the sum of n estimates, divided by n.
The Geometric Mean is the n'th root of the product of n estimates.

The average of dB power is therefore simply the Geometric Mean of the estimates, still expressed in dB.

Maybe that is not what you want to do. It is certainly not traditional power spectrum accumulation, but in signal processing you do what works and is computationally efficient.
In this case the log computation for dB can be avoided until after the accumulation, when it is best applied for the human display.