サンプリング レート 16e3 の信号があり、その周波数範囲は 125 ~ 1000 Hz です。したがって、スペクグラムをプロットすると、すべての未使用の周波数のためにかなり小さな色範囲が得られます。
斧の制限を設定して修正しようとしましたが、うまくいきません。
未使用の周波数をカットする方法や NaN に置き換える方法はありますか?
125 Hz 未満の使用されていない周波数がまだいくつかあるため、データを 2e3 にリサンプリングしても機能しません。
ご協力いただきありがとうございます。
サンプリング レート 16e3 の信号があり、その周波数範囲は 125 ~ 1000 Hz です。したがって、スペクグラムをプロットすると、すべての未使用の周波数のためにかなり小さな色範囲が得られます。
斧の制限を設定して修正しようとしましたが、うまくいきません。
未使用の周波数をカットする方法や NaN に置き換える方法はありますか?
125 Hz 未満の使用されていない周波数がまだいくつかあるため、データを 2e3 にリサンプリングしても機能しません。
ご協力いただきありがとうございます。
specgram() がすべての作業を行っています。スペックグラム関数でaxes.pyを見ると、それがどのように機能するかがわかります。元の機能はPython27\Lib\site-packages\matplotlib\axes.py私のコンピューターにあります。
<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
これをモデルにした独自の関数を作成し、必要に応じて Pxx データをクリップする必要がある場合があります。
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 は、(len(freqs),len(bins) の形状の 2 次元配列です。
>>> Pxx.shape
(129, 311)
>>> freqs.shape
(129,)
>>> bins.shape
(311,)
>>>
これにより、Pxx と freqs が制限されます
Pxx = Pxx[(freqs >= 125) & (freqs <= 1000)]
freqs = freqs[(freqs >= 125) & (freqs <= 1000)]
これが完全な解決策です - my_specgram() - gallery の specgram_demo で使用されます。
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()
specgram() は (Pxx, freqs, bins, im) を返します。ここで、im は AxesImage インスタンス [1] です。それを使用して、プロットの制限を設定できます。
Pxx, freqs, bins, im = plt.specgram(signal, Fs=fs)
im.set_ylim((125,1000))
[1] http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.specgram
これは、プロットされる周波数の範囲を変更するhttp://matplotlib.org/examples/pylab_examples/specgram_demo.htmlの適応バージョンです。
#!/usr/bin/env python
#### from the example
####
from pylab import *
dt = 0.0005
t = arange(0.0, 20.0, dt)
s1 = sin(2*pi*100*t)
s2 = 2*sin(2*pi*400*t)
mask = where(logical_and(t>10, t<12), 1.0, 0.0)
s2 = s2 * mask
nse = 0.01*randn(len(t))
x = s1 + s2 + nse # the signal
NFFT = 1024 # the length of the windowing segments
Fs = int(1.0/dt) # the sampling frequency
ax1 = subplot(211)
plot(t, x)
subplot(212, sharex=ax1)
Pxx, freqs, bins, im = specgram(x, NFFT=NFFT, Fs=Fs, noverlap=900,
cmap=cm.gist_heat)
#### edited from the example
####
# here we get access to the axes
x1,x2,y1,y2 = axis()
# leave x range the same, change y (frequency) range
axis((x1,x2,25,500))
show()