Just a quick update.

New Python code to display spectrum also with a heat-map.

from rtlsdr import RtlSdr
from matplotlib import mlab as mlab
import matplotlib.pyplot as plt
import os
from tkinter import *


# SETTINGS
GAIN = 4
SAMPLE_RATE = 2.4e6
F_CORRECTION = 60
DATA_RES = 128
try:
    sdr = RtlSdr()
except:
    print("Error: cannot initialize the SDR device")
    exit(0)

#ASK FOR FREQUENCIES
os.system("clear")
print("Insert initial frequency [MHz]: ", end="")
CENTER_FREQ = int(input())*1e6
print("Insert final frequency [MHz]: ", end="")
END_FREQ = int(input())*1e6

#CONFIGURATION
sdr.sample_rate = SAMPLE_RATE
sdr.center_freq = CENTER_FREQ
sdr.freq_correction = F_CORRECTION
sdr.gain = GAIN
garbage = sdr.read_samples(DATA_RES * 1024)
step = 2.5e5
half_step = step/2
center_f_arr = []
n = int((END_FREQ - CENTER_FREQ)/step)


def getdata():
    pow_arr = []
    freq_arr = []
    for k in range(0,n):
        sdr.center_freq = CENTER_FREQ + (step*k)
        center_f_arr.append((CENTER_FREQ + (step*k))/1e6) #save the central frequency for later
        samples = sdr.read_samples(DATA_RES * 1024)

        power, psd_freq = mlab.psd(samples, NFFT=1024, Fs=sdr.sample_rate/1e6)
        psd_freq = psd_freq + sdr.center_freq/1e6

        for i in range(0, len(psd_freq)):
            if psd_freq[i] < (sdr.center_freq+half_step)/1e6 and psd_freq[i] > (sdr.center_freq-half_step)/1e6:
                freq_arr.append(psd_freq[i])
                pow_arr.append(power[i])

        if(len(freq_arr) > 536870000):
            print("Fatal error: trying to store too much data")
            exit(0)


    i = 0
    for e in freq_arr:
        if e in center_f_arr:
            pow_arr[i] = pow_arr[i-1]
            pow_arr[i+1] = pow_arr[i+2]
            pow_arr[i-1] = pow_arr[i-2]
        i=i+1
    return freq_arr,pow_arr
 

#ask for data
freq_arr,pow_arr = getdata()

#create 2d map
mappa = [[0 for x in range(len(pow_arr))] for y in range(len(pow_arr))] 
for i in range(0, len(pow_arr)):
    for j in range(0, len(pow_arr)):
        mappa[i][j] = pow_arr[j]


f, axarr = plt.subplots(2)
axarr[0].plot(freq_arr, pow_arr, color="green")
c = axarr[1].imshow(mappa, cmap='hot', interpolation='nearest',origin ='lower', aspect='auto')
axarr[1].axis("off")
f.colorbar(c)
plt.show()

Output plot: