How To Plot Fft Of Signal With Correct Frequencies On X-axis?

I can plot signals I receive from a RTL-SDR with Matplotlib's plt.psd() method, which results in the following plot: The ultimate goal of what I'm trying to achieve is to retrieve

Solution 1:

I think your problem comes from taking the FFT of the whole signal which results in a too high frequency resolution that causes the noise you see. Matplotlib psd breaks the signal in shorter overlapping blocks, calculates the FFT of each block and averages. The function welch in Scipy signal also does this. You will get a spectrum centered around 0 Hz. You can then offset the returned frequency vector to get your original frequency range by adding the center frequency to the frequency vector.

When you use welch, the returned frequency and power vectors are not sorted in ascending frequency order. You must fftshift the output before you plot. Also, the sample frequency you pass welch must be a float. Make sure to use the scaling="spectrum" option to get the power instead of the power density. To get the correct power value, you also need to scale the power to account for the windowing effect. For a hann window, you need to divide by 1.5. Here is a working example:

from scipy.signal import welch
from numpy.fft import fftshift, fft
from numpy import arange, exp, pi, mean
import matplotlib.pyplot as plt

#generate a 1000 Hz complex tone singnal
sample_rate = 48000.#sample rate must be a float
t=arange(1024*1024)/sample_rate
signal=exp(1j*2000*pi*t)
#amplitude correction factor
corr=1.5#calculate the psd with welch
sample_freq, power = welch(signal, fs=sample_rate, window="hann", nperseg=256, noverlap=128, scaling='spectrum')

#fftshift the output 
sample_freq=fftshift(sample_freq)
power=fftshift(power)/corr

#check that the power sum is rightprintsum(power)

plt.figure(figsize=(9.84, 3.94))
plt.plot(sample_freq, power)
plt.xlabel("Frequency (MHz)")
plt.ylabel("Relative power (dB)")  
plt.show()

edit

I see three reasons that can explain why you don't get the same amplitude as the Matplotlib PSD function. By the way if you look at the doc, Matplotlib PSD has three return arguments : PSD, frequency vector and the line object so if you want the same PSD as you get by the Matplotlib function, you can just grab the data there. But I suggest you to read further to make sure you know what you're doing (Just because Matplotlib returns a value to you doesn't mean it's right or it's the value you need).

1. The way you calculate decibels is wrong. First, plotting in decibels is not the same as plotting on a logarithmic axis. The decibel scale is logarithmic while the data you plot on the logarithmic graph is still linear so obviously you won't get the same values. The decibel scale is relative which means that you compare your value to a reference. The way to calculate decibels depend on the type of units you're dealing with. In the case of primary physical units (Volts, Pascals, m/s, etc.), you would do, if we assume the units are volts: 20*log10(V/Vref) where Vref is the reference value. Now if you're dealing with squared units: energy, intensity, power density, etc. you do: 10*log10(P/Pref) where P is the squared quantity. This is because squaring in the linear domain is equivalent to a multiplication by two in the log domain. Anyway in your case the 10*log10 form applies. As for the refrence value, it is arbitrary and is often a convention in a particular field. For example, the international refrence for sound pressure in acoustics is 2e-5 Pa. Since the reference is arbitrary, you can just set it to 1 when you don't have any standard value to compare against.

2. Second, the default for Matplotlib psd is to have the "scale_by_freq" option set to true. This gives you the power by Hz or in your case by Mhz. On the other hand, the spectrum option in welch gives you the power per band. So in Matplotlib, the power is divided by the frequecy range (2.8 MHz) while in welch it's divided by the number of bands (2048). If I take the decibel ratio of the two I get 10*log10(2048/2.8)=28.6 dB which seems quite close to the difference you get.

3. Finally the correction factor you use is different depending on what you want to achieve. The reason you need a correction factor in the first place is because of the modification to the signal introduced by the windowing function. The window effectively reduces the total energy of the signal. You the need to multiply by a correction to get the right energy back. The windowing also affects the spectrum by spreading the energy to the adjacent bands. The consequence of this is that the height of the peaks in the signal is modified. Therefore, if you want the right peak height, you have to use a correction factor, while if you want the correct energy, you use another. The ratio between the two for hann window is 1.5. Usually, the amplitude correction (correct peak height) is used for display purposes while the energy correction is used when total energy is important. I would think Matplotlib PSD is amplitude corrected while the example I gave you is energy corrected so there could be a 1.5 factor there.

Since what you want to do is find peaks in the spectrum, maybe the whole amplitude thing is not so important after all. Peak finding usually works with the relative difference between frequency bands.