[SciPy-Dev] 2D histogram: request for plotting variable bin size

[Sturla Molden]

Hi Frank

I use Python to plot wavelet spectrograms with logarithmically scaled 
frequencies. I e.g. interpolate the spectrum to the frequencies I like 
to plot, and pass that to pcolor or contourf. Usually I plot with a 
linear y-axis first, and then use set_yscale('log'). Usually the y-axis 
must be flipped too.

Then I e.g. end up with something like this, plotting with 250 
logartihmically scaled frequencies from 10 Hz to 10 kHz:

cwt,scale= wavelet(signal, sampling_rate)
time = np.arange(signal.shape[0])/sampling_rate
power = np.abs(cwt)**2/scale # bias corrected power
periods = wavelet_to_fourier(scale)
spectrum_freqs = np.logspace(np.log10(10), np.log10(10000), 250)
spectrum_periods = spectrum_freqs**-1
spectrum = interpolate_spectrum(spectrum_periods, power, periods)
X,Y = np.meshgrid(time, spectrum_freqs)
cs = plt.contourf(X,Y,spectrum, cmap=cm.jet, origin=origin)
ax = plt.gca()
ax.set_yscale('log')
ylim = ax.get_ylim()
ax.set_ylim(ylim[::-1])



import numpy as np
from bisect import bisect_left

def interpolate_spectrum(spectrum_periods, wavelet_power, wavelet_periods):

     """ linear interpolation for plotting spectrum at given frequencies 
"""
     n0 = spectrum_periods.shape[0]
     n1 = wavelet_power.shape[1]
     m = wavelet_power.shape[0]
     spectrum_power = np.zeros((n0,n1))
     for k,sp in enumerate(spectrum_periods):
         i = bisect_left(wavelet_periods,sp)
         j = i-1
         i = 0 if i<0 else i
         i = m-1 if i > m-1 else i
         j = 0 if j<0 else j
         j = m-1 if j > m-1 else j
         eps = 1E-5 # avoid log of zero
         di = 1./(np.abs(sp - wavelet_periods[i]) + eps)
         dj = 1./(np.abs(sp - wavelet_periods[j]) + eps)
         spectrum_power[k,:] = (di*wavelet_power[i,:] + 
dj*wavelet_power[j,:])/(di+dj)
     return spectrum_power



Not sure if this helps or not, but it might.


Sturla



On 31.01.2013 16:00, Frank Breitling wrote:
> Hi Skipper,
>
> I would like to use Python for plotting dynamic radio spectra e.g. as
> found at
>
> http://science.nasa.gov/media/medialibrary/1998/12/20/ast22dec98_1_resources/leonid_dynspectrum.jpg
>
> (Image is also attached).
> Given the high resolution in time and frequency I believe broken_barh is
> not suited.
>
> But thanks for your suggestion
>
> Frank
>
>
> On 31.01.2013 15:17, Skipper Seabold wrote:
>> On Thu, Jan 31, 2013 at 9:01 AM, Frank Breitling <fbreitling at aip.de
>> <mailto:fbreitling at aip.de>> wrote:
>>
>>     Hello,
>>
>>     http://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram2d.html
>>
>>     explains the usage of 2D histograms.
>>     Although it is possible to construct 2D histograms with a variable bin
>>     size, there is no simple way to plot them.
>>     Therefore I would like to request an implementation to imshow that
>>     allows the visualization of 2D histograms with variable bin size. I
>>     imagine a syntax like
>>
>>     imshow(historgram2d, xedges=xedges, yedges=yedges)
>>
>>     where xedges and yedges are the bin edges along the x- and y-axis.
>>
>>     I have attached a short program and its output which constructs a
>>     numpy.historgram2d with variable bin width (xedges=[0,1,3]) and shows
>>     its bin content with imshow.
>>     I would highly appreciate if imshow would be extended to represent the
>>     histogram correctly.
>>     If there is an alternative solution, I would be interested to learn
>>     about it as well.
>>
>>
>> This may be a discussion for matplotlib-user, but have you had a look
>> at trying to use broken_barh for your needs?
>>
>> http://matplotlib.org/examples/pylab_examples/broken_barh.html
>>
>> Skipper
>>
>>
>> _______________________________________________
>> SciPy-Dev mailing list
>> SciPy-Dev at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>