Cutting of unused frequencies in specgram matplotlib

Question

I have a signal with sampling rate 16e3, its frequency range is from 125 to 1000 Hz. So if i plot a specgram i get a pretty small colorrange because of all the unused frequencys.

ive tried to fix it with setting ax limits but that does not work.

is there any way to cut off unused frequencys or replace them with NaNs?

Resampling the Data to 2e3 won't work because there are still some not used frequencys below 125 Hz.

Thanks for your help.

Solution 1

specgram() is doing all the work for you. If you look in axes.py at the specgram function you can see how it works. The original function is in Python27\Lib\site-packages\matplotlib\axes.py on my computer.

    <snip>  
    Pxx, freqs, bins = mlab.specgram(x, NFFT, Fs, detrend,
         window, noverlap, pad_to, sides, scale_by_freq)

    Z = 10. * np.log10(Pxx)
    Z = np.flipud(Z)

    if xextent is None: xextent = 0, np.amax(bins)
    xmin, xmax = xextent
    freqs += Fc
    extent = xmin, xmax, freqs[0], freqs[-1]
    im = self.imshow(Z, cmap, extent=extent, **kwargs)
    self.axis('auto')

    return Pxx, freqs, bins, im

You might have to make your own function modeled on this and clip the Pxx data to suit your needs.

    Pxx, freqs, bins = mlab.specgram(x, NFFT, Fs, detrend,
         window, noverlap, pad_to, sides, scale_by_freq)
    # ****************
    # create a new limited Pxx and freqs
    # 
    # ****************
    Z = 10. * np.log10(Pxx)
    Z = np.flipud(Z)

Pxx is a 2d array with a shape of (len(freqs),len(bins)

>>> Pxx.shape
(129, 311)
>>> freqs.shape
(129,)
>>> bins.shape
(311,)
>>> 

This will limit Pxx and freqs

Pxx = Pxx[(freqs >= 125) & (freqs <= 1000)]
freqs = freqs[(freqs >= 125) & (freqs <= 1000)]

Here is a complete solution - my_specgram() - used with the specgram_demo from the gallery.

from pylab import *
from matplotlib import *


# 100, 400 and 200 Hz sine 'wave'
dt = 0.0005
t = arange(0.0, 20.0, dt)
s1 = sin(2*pi*100*t)
s2 = 2*sin(2*pi*400*t)
s3 = 2*sin(2*pi*200*t)

# create a transient "chirp"
mask = where(logical_and(t>10, t<12), 1.0, 0.0)
s2 = s2 * mask

# add some noise into the mix
nse = 0.01*randn(len(t))

x = s1 + s2 + +s3 + nse # the signal
NFFT = 1024       # the length of the windowing segments
Fs = int(1.0/dt)  # the sampling frequency

# modified specgram()
def my_specgram(x, NFFT=256, Fs=2, Fc=0, detrend=mlab.detrend_none,
             window=mlab.window_hanning, noverlap=128,
             cmap=None, xextent=None, pad_to=None, sides='default',
             scale_by_freq=None, minfreq = None, maxfreq = None, **kwargs):
    """
    call signature::

      specgram(x, NFFT=256, Fs=2, Fc=0, detrend=mlab.detrend_none,
               window=mlab.window_hanning, noverlap=128,
               cmap=None, xextent=None, pad_to=None, sides='default',
               scale_by_freq=None, minfreq = None, maxfreq = None, **kwargs)

    Compute a spectrogram of data in *x*.  Data are split into
    *NFFT* length segments and the PSD of each section is
    computed.  The windowing function *window* is applied to each
    segment, and the amount of overlap of each segment is
    specified with *noverlap*.

    %(PSD)s

      *Fc*: integer
        The center frequency of *x* (defaults to 0), which offsets
        the y extents of the plot to reflect the frequency range used
        when a signal is acquired and then filtered and downsampled to
        baseband.

      *cmap*:
        A :class:`matplotlib.cm.Colormap` instance; if *None* use
        default determined by rc

      *xextent*:
        The image extent along the x-axis. xextent = (xmin,xmax)
        The default is (0,max(bins)), where bins is the return
        value from :func:`mlab.specgram`

      *minfreq, maxfreq*
        Limits y-axis. Both required

      *kwargs*:

        Additional kwargs are passed on to imshow which makes the
        specgram image

      Return value is (*Pxx*, *freqs*, *bins*, *im*):

      - *bins* are the time points the spectrogram is calculated over
      - *freqs* is an array of frequencies
      - *Pxx* is a len(times) x len(freqs) array of power
      - *im* is a :class:`matplotlib.image.AxesImage` instance

    Note: If *x* is real (i.e. non-complex), only the positive
    spectrum is shown.  If *x* is complex, both positive and
    negative parts of the spectrum are shown.  This can be
    overridden using the *sides* keyword argument.

    **Example:**

    .. plot:: mpl_examples/pylab_examples/specgram_demo.py

    """

    #####################################
    # modified  axes.specgram() to limit
    # the frequencies plotted
    #####################################

    # this will fail if there isn't a current axis in the global scope
    ax = gca()
    Pxx, freqs, bins = mlab.specgram(x, NFFT, Fs, detrend,
         window, noverlap, pad_to, sides, scale_by_freq)

    # modified here
    #####################################
    if minfreq is not None and maxfreq is not None:
        Pxx = Pxx[(freqs >= minfreq) & (freqs <= maxfreq)]
        freqs = freqs[(freqs >= minfreq) & (freqs <= maxfreq)]
    #####################################

    Z = 10. * np.log10(Pxx)
    Z = np.flipud(Z)

    if xextent is None: xextent = 0, np.amax(bins)
    xmin, xmax = xextent
    freqs += Fc
    extent = xmin, xmax, freqs[0], freqs[-1]
    im = ax.imshow(Z, cmap, extent=extent, **kwargs)
    ax.axis('auto')

    return Pxx, freqs, bins, im

# plot
ax1 = subplot(211)
plot(t, x)
subplot(212, sharex=ax1)

# the minfreq and maxfreq args will limit the frequencies 
Pxx, freqs, bins, im = my_specgram(x, NFFT=NFFT, Fs=Fs, noverlap=900, 
                                cmap=cm.Accent, minfreq = 180, maxfreq = 220)
show()
close()

Solution 2

These days, there's an easier way to do this than when the question was asked: you can use matplotlib.pyplot.axis to set ymin and ymax to your desired frequencies. It's quite easy; here's a snippet from my program:

plt.specgram(xmit, NFFT=65536, Fs=Fs)
plt.axis(ymin=Fc-Fa*10, ymax=Fc+Fa*10)
plt.show()

Comments

[Admin]

thanks alot. would upvote if i had enough reputation.

[wwii]

I added my_specgram() which I think does what you need.

[user7610]

now you should have enough reputation to give that upvote ;-)

[moeseth]

@wwii, hi, why do you do Z = 10. * np.log10(Pxx) to spectrum?

[wwii]

@moeseth That is part of the original function - I bracketed my changes, above and below, with a line of hashes.

[wwii]

@moeseth one reason for using logs is to compress information to a more usable (viewable(?)) form.

[wwii]

@moeseth - looks like the docstring my current version of axes.py has this addition - Also note that while the plot is in dB, the *Pxx* array returned is linear in power. So they are compressing the data and displaying it as dB power which is common for spectral "graphs".