def pure_data_plot(self,connect=False,suffix='',cmap=cm.jet,bg=cm.bone(0.3)):
     #fig=plt.figure()
     ax=plt.axes()
     plt.axhline(y=0,color='grey', zorder=-1)
     plt.axvline(x=0,color='grey', zorder=-2)
     if cmap is None:
         if connect: ax.plot(self.x,self.y, 'b-',lw=2,alpha=0.5)
         ax.scatter(self.x,self.y, marker='o', c='b', s=40)
     else:
         if connect:
             if cmap in [cm.jet,cm.brg]:
                 ax.plot(self.x,self.y, 'c-',lw=2,alpha=0.5,zorder=-1)
             else:
                 ax.plot(self.x,self.y, 'b-',lw=2,alpha=0.5)
         c=[cmap((f-self.f[0])/(self.f[-1]-self.f[0])) for f in self.f]
         #c=self.f
         ax.scatter(self.x, self.y, marker='o', c=c, edgecolors=c, zorder=True, s=40) #, cmap=cmap)
     #plt.axis('equal')
     ax.set_xlim(xmin=-0.2*amax(self.x), xmax=1.2*amax(self.x))
     ax.set_aspect('equal')  #, 'datalim')
     if cmap in [cm.jet,cm.brg]:
         ax.set_axis_bgcolor(bg)
     if self.ZorY == 'Z':
         plt.xlabel(r'resistance $R$ in Ohm'); plt.ylabel(r'reactance $X$ in Ohm')
     if self.ZorY == 'Y':
         plt.xlabel(r'conductance $G$ in Siemens'); plt.ylabel(r'susceptance $B$ in Siemens')
     if self.show: plt.show()
     else: plt.savefig(join(self.sdc.plotpath,'c{}_{}_circle_data'.format(self.sdc.case,self.ZorY)+self.sdc.suffix+self.sdc.outsuffix+suffix+'.png'), dpi=240)
     plt.close()

def plot_zipf(*freq):
	'''
	basic plotting using matplotlib and pylab
	'''
	ranks, frequencies = [], []
	langs, colors = [], []
	langs = ["English", "German", "Finnish"]
	colors = ['#FF0000', '#00FF00', '#0000FF']
	if bonus_part:
		colors.extend(['#00FFFF', '#FF00FF', '#FFFF00'])
		langs.extend(["English (Stemmed)", "German (Stemmed)", "Finnish (Stemmed)"])

	plt.subplot(111) # 1, 1, 1

	num = 6 if bonus_part else 3
	for i in xrange(num):
		ranks.append(range(1, len(freq[i]) + 1))
		frequencies.append([e[1] for e in freq[i]])

		# log x and y axi, both with base 10
		plt.loglog(ranks[i], frequencies[i], marker='', basex=10, color=colors[i], label=langs[i])

	plt.legend()
	plt.grid(True)
	plt.title("Zipf's law!")

	plt.xlabel('Rank')
	plt.ylabel('Frequency')

	plt.show()

    def draw(cls, t_max, agents_proportions, eco_idx, parameters):

        color_set = ["green", "blue", "red"]

        for agent_type in range(3):
            plt.plot(np.arange(t_max), agents_proportions[:, agent_type],
                     color=color_set[agent_type], linewidth=2.0, label="Type-{} agents".format(agent_type))

            plt.ylim([-0.1, 1.1])

        plt.xlabel("$t$")
        plt.ylabel("Proportion of indirect exchanges")

        # plt.suptitle('Direct choices proportion per type of agents', fontsize=14, fontweight='bold')
        plt.legend(loc='upper right', fontsize=12)

        print(parameters)

        plt.title(
            "Workforce: {}, {}, {};   displacement area: {};   vision area: {};   alpha: {};   tau: {}\n"
            .format(
                parameters["x0"],
                parameters["x1"],
                parameters["x2"],
                parameters["movement_area"],
                parameters["vision_area"],
                parameters["alpha"],
                parameters["tau"]
                          ), fontsize=12)

        if not path.exists("../../figures"):
            mkdir("../../figures")

        plt.savefig("../../figures/figure_{}.pdf".format(eco_idx))
        plt.show()

def example_filterbank():
    from pylab import plt
    import numpy as np

    x = _create_impulse(2000)
    gfb = GammatoneFilterbank(density=1)

    analyse = gfb.analyze(x)
    imax, slopes = gfb.estimate_max_indices_and_slopes()
    fig, axs = plt.subplots(len(gfb.centerfrequencies), 1)
    for (band, state), imx, ax in zip(analyse, imax, axs):
        ax.plot(np.real(band))
        ax.plot(np.imag(band))
        ax.plot(np.abs(band))
        ax.plot(imx, 0, 'o')
        ax.set_yticklabels([])
        [ax.set_xticklabels([]) for ax in axs[:-1]]

    axs[0].set_title('Impulse responses of gammatone bands')

    fig, ax = plt.subplots()

    def plotfun(x, y):
        ax.semilogx(x, 20*np.log10(np.abs(y)**2))

    gfb.freqz(nfft=2*4096, plotfun=plotfun)
    plt.grid(True)
    plt.title('Absolute spectra of gammatone bands.')
    plt.xlabel('Normalized Frequency (log)')
    plt.ylabel('Attenuation /dB(FS)')
    plt.axis('Tight')
    plt.ylim([-90, 1])
    plt.show()

    return gfb

def serve_css(name, length, keys, values):
    from pylab import plt, mpl
    mpl.rcParams['font.sans-serif'] = ['SimHei']
    mpl.rcParams['axes.unicode_minus'] = False
    from matplotlib.font_manager import FontProperties
    # font = FontProperties(fname="d:\Users\ll.tong\Desktop\msyh.ttf", size=12)
    font = FontProperties(fname="/usr/share/fonts/msyh.ttf", size=11)
    plt.xlabel(u'')
    plt.ylabel(u'出现次数',fontproperties=font)
    plt.title(u'词频统计',fontproperties=font)
    plt.grid()
    keys = keys.decode("utf-8").split(' ')
    values = values.split(' ')
    valuesInt = []
    for value in values:
        valuesInt.append(int(value))

    plt.xticks(range(int(length)), keys)
    plt.plot(range(int(length)), valuesInt)
    plt.xticks(rotation=defaultrotation, fontsize=9,fontproperties=font)
    plt.yticks(fontsize=10,fontproperties=font)
    name = name + str(datetime.now().date()).replace(':', '') + '.png'
    imgUrl = 'static/temp/' + name
    fig = matplotlib.pyplot.gcf()
    fig.set_size_inches(12.2, 2)
    plt.savefig(imgUrl, bbox_inches='tight', figsize=(20,4), dpi=100)
    plt.close()
    tempfile = static_file(name, root='./static/temp/')
    #os.remove(imgUrl)
    return tempfile

def plot_weightings():
    """Plots all weighting functions defined in :module: splweighting."""
    from scipy.signal import freqz
    from pylab import plt, np

    sample_rate = 48000
    num_samples = 2*4096

    fig, ax = plt.subplots()

    for name, weight_design in sorted(
            _weighting_coeff_design_funsd.items()):
        b, a = weight_design(sample_rate)
        w, H = freqz(b, a, worN=num_samples)

        freq = w*sample_rate / (2*np.pi)

        ax.semilogx(freq, 20*np.log10(np.abs(H)+1e-20),
                    label='{}-Weighting'.format(name))

    plt.legend(loc='lower right')
    plt.xlabel('Frequency / Hz')
    plt.ylabel('Damping / dB')
    plt.grid(True)
    plt.axis([10, 20000, -80, 5])
    return fig, ax

def plot_response(data, plate_name, save_folder = 'Figures/'):
    """
    """
    if not os.path.isdir(save_folder):
        os.makedirs(save_folder)

    for block in data:
        #
        group = group_similar(data[block].keys())
        names = data[block].keys()
        names.sort()
        #
        plt.figure(figsize=(16, 4 + len(names)/8), dpi=300)
        #
        for i, name in enumerate(names):
            a, b, c = get_index(group, name)
            color, pattern = color_shade_pattern(a, b, c, group)
            mean = data[block][name]['mean'][0]
            std = data[block][name]['std'][0]

            plt.barh([i], [mean], height=1.0, color=color, hatch=pattern)
            plt.errorbar([mean], [i+0.5], xerr=[std], ecolor = [0,0,0], linestyle = '')

        plt.yticks([i+0.5 for i in xrange(len(names))], names, size = 8)
        plt.title(plate_name)
        plt.ylim(0, len(names))
        plt.xlabel('change')
        plt.tight_layout()

        plt.savefig(save_folder + 'response_' + str(block + 1))
    #
    return None

    def generate_start_time_figures(self):
        recording_time_grouped_by_patient = self.pain_data[["PatientID", "NRSTimeFromEndSurgery_mins"]].groupby("PatientID")
        recording_start_minutes = recording_time_grouped_by_patient.min()

        fig1 = "fig1.pdf"
        fig2 = "fig2.pdf"

        plt.figure(figsize=[8,4])
        plt.title("Pain score recording start times", fontsize=14).set_y(1.05) 
        plt.ylabel("Occurrences", fontsize=14)
        plt.xlabel("Recording Start Time (minutes)", fontsize=14)
        plt.hist(recording_start_minutes.values, bins=20, color="0.5")
        plt.savefig(os.path.join(self.tmp_directory, fig1), bbox_inches="tight")

        plt.figure(figsize=[8,4])
        plt.title("Pain score recording start times, log scale", fontsize=14).set_y(1.05) 
        plt.ylabel("Occurrences", fontsize=14)
        plt.xlabel("Recording Start Time (minutes)", fontsize=14)
        plt.hist(recording_start_minutes.values, bins=20, log=True, color="0.5")
        plt.savefig(os.path.join(self.tmp_directory, fig2), bbox_inches="tight")

        #save the figures in panel format
        f = open(os.path.join(self.tmp_directory, "tmp.tex"), 'w')
        f.write(r"""
            \documentclass[%
            ,float=false % this is the new default and can be left away.
            ,preview=true
            ,class=scrartcl
            ,fontsize=20pt
            ]{standalone}
            \usepackage[active,tightpage]{preview}
            \usepackage{varwidth}
            \usepackage{graphicx}
            \usepackage[justification=centering]{caption}
            \usepackage{subcaption}
            \usepackage[caption=false,font=footnotesize]{subfig}
            \renewcommand{\thesubfigure}{\Alph{subfigure}}
            \begin{document}
            \begin{preview}
            \begin{figure}[h]
                \begin{subfigure}{0.5\textwidth}
                        \includegraphics[width=\textwidth]{""" + fig1 + r"""}
                        \caption{Normal scale}
                \end{subfigure}\begin{subfigure}{0.5\textwidth}
                        \includegraphics[width=\textwidth]{""" + fig2 + r"""}
                        \caption{Log scale}
                \end{subfigure}
            \end{figure}
            \end{preview}
            \end{document}
        """)
        f.close()
        subprocess.call(["pdflatex", 
                            "-halt-on-error", 
                            "-output-directory", 
                            self.tmp_directory, 
                            os.path.join(self.tmp_directory, "tmp.tex")])
        shutil.move(os.path.join(self.tmp_directory, "tmp.pdf"), 
                    os.path.join(self.output_directory, "pain_score_start_times.pdf"))

    def plot(self, new_plot=False, xlim=None, ylim=None, title=None, figsize=None,
             xlabel=None, ylabel=None, fontsize=None, show_legend=True, grid=True):
        """
        Plot data using matplotlib library. Use show() method for matplotlib to see result or ::

            %pylab inline

        in IPython to see plot as cell output.

        :param bool new_plot: create or not new figure
        :param xlim: x-axis range
        :param ylim: y-axis range
        :type xlim: None or tuple(x_min, x_max)
        :type ylim: None or tuple(y_min, y_max)
        :param title: title
        :type title: None or str
        :param figsize: figure size
        :type figsize: None or tuple(weight, height)
        :param xlabel: x-axis name
        :type xlabel: None or str
        :param ylabel: y-axis name
        :type ylabel: None or str
        :param fontsize: font size
        :type fontsize: None or int
        :param bool show_legend: show or not labels for plots
        :param bool grid: show grid or not

        """
        xlabel = self.xlabel if xlabel is None else xlabel
        ylabel = self.ylabel if ylabel is None else ylabel
        figsize = self.figsize if figsize is None else figsize
        fontsize = self.fontsize if fontsize is None else fontsize
        self.fontsize_ = fontsize
        self.show_legend_ = show_legend
        title = self.title if title is None else title
        xlim = self.xlim if xlim is None else xlim
        ylim = self.ylim if ylim is None else ylim
        new_plot = self.new_plot or new_plot

        if new_plot:
            plt.figure(figsize=figsize)

        plt.xlabel(xlabel, fontsize=fontsize)
        plt.ylabel(ylabel, fontsize=fontsize)
        plt.title(title, fontsize=fontsize)
        plt.tick_params(axis='both', labelsize=fontsize)
        plt.grid(grid)

        if xlim is not None:
            plt.xlim(xlim)

        if ylim is not None:
            plt.ylim(ylim)

        self._plot()

        if show_legend:
            plt.legend(loc='best', scatterpoints=1)

 def plot_smoothed_alpha_comparison(self,rmsval,suffix=''):
     plt.plot(self.f,self.alpha,'ko',label='data set')
     plt.plot(self.f,self.salpha,'c-',lw=2,label='smoothed angle $\phi$')
     plt.xlabel('frequency in Hz')
     plt.ylabel('angle $\phi$ in coordinates of circle')
     plt.legend()
     ylims=plt.axes().get_ylim()
     plt.yticks((arange(9)-4)*0.5*pi, ['$-2\pi$','$-3\pi/2$','$-\pi$','$-\pi/2$','$0$','$\pi/2$','$\pi$','$3\pi/2$','$2\pi$'])
     plt.ylim(ylims)
     plt.title('RMS offset from smooth curve: {:.4f}'.format(rmsval))
     if self.show: plt.show()
     else: plt.savefig(join(self.sdc.plotpath,'salpha','c{}_salpha_on_{}_circle'.format(self.sdc.case,self.ZorY)+self.sdc.suffix+self.sdc.outsuffix+suffix+'.png'), dpi=240)
     plt.close()

def convert_all_to_png(vis_path, out_dir = "maps_png", size = None) :

    units = { 'gas_density' : 'Gas Density [g/cm$^3$]',
              'Tm' : 'Temperature [K]',
              'Tew' : 'Temperature [K]',
              'S' : 'Entropy []',
              'dm' : 'DM Density [g/cm$^3$]',
              'v' : 'Velocity [km/s]' }

    log_list = ['gas_density']

    for vis_file in os.listdir(vis_path) :
        if ".dat" not in vis_file :
            continue
        print "converting %s" % vis_file
        map_type = re.search('sigma_(.*)_[xyz]', vis_file).group(1)

        (image, pixel_size, axis_values) = read_visualization_data(vis_path+"/"+vis_file, size)
        print "image width in Mpc/h: ", axis_values[-1]*2.0

        x, y = np.meshgrid( axis_values, axis_values )

        cmap_max = image.max()
        cmap_min = image.min()


        ''' plotting '''
        plt.figure(figsize=(5,4))

        if map_type in log_list:
            plt.pcolor(x,y,image, norm=LogNorm(vmax=cmap_max, vmin=cmap_min))
        else :
            plt.pcolor(x,y,image, vmax=cmap_max, vmin=cmap_min)

        cbar = plt.colorbar()
        if map_type in units.keys() :
            cbar.ax.set_ylabel(units[map_type])

        plt.axis([axis_values[0], axis_values[-1],axis_values[0], axis_values[-1]])

        del image

        plt.xlabel(r"$Mpc/h$", fontsize=18)
        plt.ylabel(r"$Mpc/h$", fontsize=18)

        out_file = vis_file.replace("dat", "png")

        plt.savefig(out_dir+"/"+out_file, dpi=150 )

        plt.close()
        plt.clf()

def plot_fault_framework(fault_framework):
    fm = fault_framework
    plt.plot(fm.Y_PC, fm.DEP, '-o')
    plt.axis('equal')
    plt.axhline(0, color='black')
    plt.gca().set_yticks(fm.DEP)
    plt.gca().set_xticks(fm.Y_PC)
    plt.grid('on')
    plt.xlabel('From trench to continent(km)')
    plt.ylabel('depth (km)')

    for xi, yi, dip in zip(fm.Y_PC, fm.DEP, fm.DIP_D):
        plt.text(xi, yi, 'dip = %.1f'%dip)

    plt.gca().invert_yaxis()

 def plot_overview(self,suffix=''):
     x=self.x; y=self.y; r=self.radius; cx,cy=self.center.real,self.center.imag
     ax=plt.axes()
     plt.scatter(x,y, marker='o', c='b', s=40)
     plt.axhline(y=0,color='grey', zorder=-1)
     plt.axvline(x=0,color='grey', zorder=-2)
     t=linspace(0,2*pi,201)
     circx=r*cos(t) + cx
     circy=r*sin(t) + cy
     plt.plot(circx,circy,'g-')
     plt.plot([cx],[cy],'gx',ms=12)
     if self.ZorY == 'Z':
         philist,flist=[self.phi_a,self.phi_p,self.phi_n],[self.fa,self.fp,self.fn]
     elif self.ZorY == 'Y':
         philist,flist=[self.phi_m,self.phi_s,self.phi_r],[self.fm,self.fs,self.fr]
     for p,f in zip(philist,flist):
         if f is not None:
             xpos=cx+r*cos(p); ypos=cy+r*sin(p); xos=0.2*(xpos-cx); yos=0.2*(ypos-cy)
             plt.plot([0,xpos],[0,ypos],'co-')
             ax.annotate('{:.3f} Hz'.format(f), xy=(xpos,ypos),  xycoords='data',
                         xytext=(xpos+xos,ypos+yos), textcoords='data', #textcoords='offset points',
                         arrowprops=dict(arrowstyle="->", shrinkA=0, shrinkB=10)
                         )
     #plt.xlim(0,0.16)
     #plt.ylim(-0.1,0.1)
     plt.axis('equal')
     if self.ZorY == 'Z':
         plt.xlabel(r'resistance $R$ in Ohm'); plt.ylabel(r'reactance $X$ in Ohm')
     if self.ZorY == 'Y':
         plt.xlabel(r'conductance $G$ in Siemens'); plt.ylabel(r'susceptance $B$ in Siemens')
     plt.title("fitting the admittance circle with Powell's method")
     tx1='best fit (fmin_powell):\n'
     tx1+='center at G+iB = {:.5f} + i*{:.8f}\n'.format(cx,cy)
     tx1+='radius = {:.5f};  '.format(r)
     tx1+='residue: {:.2e}'.format(self.resid)
     txt1=plt.text(-r,cy-1.1*r,tx1,fontsize=8,ha='left',va='top')
     txt1.set_bbox(dict(facecolor='gray', alpha=0.25))
     idxlist=self.to_be_annotated('triple')
     ofs=self.annotation_offsets(idxlist,factor=0.1,xshift=0.15)
     for i,j in enumerate(idxlist):
         xpos,ypos = x[j],y[j]; xos,yos = ofs[i].real,ofs[i].imag
         ax.annotate('{:.1f} Hz'.format(self.f[j]), xy=(xpos,ypos),  xycoords='data',
                     xytext=(xpos+xos,ypos+yos), textcoords='data', #textcoords='offset points',
                     arrowprops=dict(arrowstyle="->", shrinkA=0, shrinkB=10)
                     )
     if self.show: plt.show()
     else: plt.savefig(join(self.sdc.plotpath,'c{}_fitted_{}_circle'.format(self.sdc.case,self.ZorY)+suffix+'.png'), dpi=240)
     plt.close()

    def test_dep(self):
        xf = arange(0, 425)
        deps = self.fm.get_dep(xf)
        plt.plot(xf,deps)

        plt.gca().set_yticks(self.fm.DEP)
        plt.gca().set_xticks(self.fm.Y_PC)
        
        plt.grid('on')
        plt.title('Ground x versus depth')
        plt.xlabel('Ground X (km)')
        plt.ylabel('depth (km)')
        plt.axis('equal')
        plt.gca().invert_yaxis()
        plt.savefig(join(self.outs_dir, '~Y_PC_vs_deps.png'))
        plt.close()

    def update_img((expected, output)):
        plt.cla()
        plt.ylim((vmin, vmin+vmax))
        plt.xlim((vmin, vmin+vmax))
        ax = fig.add_subplot(111)
        plt.plot([vmin, vmin+vmax], [vmin, vmin+vmax])
        ax.grid(True)
        plt.xlabel("expected output")
        plt.ylabel("network output")
        plt.legend()

        expected = expected*vmax + vmin
        output = output*vmax + vmin
        #scat.set_offsets((expected, output))
        scat = ax.scatter(expected, output)
        return scat

def plot_baseline(data, plate_name, save_folder = r'Figures/'):
    """
    """
    colors = ((0.2, 0.2, 0.2),
              (0.5, 0.5, 0.5),
              (0.7, 0.7, 0.7),
              (0.3, 0.3, 0.3))

    names = data.keys()
    names.sort()
    fig, axs = plt.subplots(figsize=(8,3))
    for index, name in enumerate(names):
        for value in data[name]['original_data']:
            plot_color = colors[index % len(colors)]
            
            if abs(value - data[name]['mean'][0]) > data[name]['std'][0] * 2.0:
                axs.plot([value], [index], 'ko', markerfacecolor = [1,1,1])
            else:
                axs.plot([value], [index], 'ko', color = plot_color)

        axs.plot([data[name]['mean'][0] for _ in xrange(2)],
                 [index-0.25, index+0.25],
                  'k-')
        axs.plot([data[name]['mean'][0] - data[name]['std'][0] for _ in xrange(2)],
                 [index-0.25, index+0.25],
                 'k--')
        axs.plot([data[name]['mean'][0] + data[name]['std'][0] for _ in xrange(2)],
                 [index-0.25, index+0.25],
                 'k--')

    plt.yticks([i for i in xrange(len(names))], names, size = 10)
    plt.title(plate_name)
    plt.ylim(-0.5,len(names)-0.5)
    plt.xlabel('Fluorescent intensity')
    plt.tight_layout()

    save_filename = save_folder + 'baseline_average'

    pdf = PdfPages(save_filename.split('.')[0] + '.pdf')
    pdf.savefig(fig)
    pdf.close()

    plt.savefig(save_filename)
    #
    return None

 def start_plot(self,w=1.3,connect=False):
     self.fig=plt.figure()
     self.ax=plt.axes()
     plt.axhline(y=0,color='grey', zorder=-1)
     plt.axvline(x=0,color='grey', zorder=-2)
     self.plot_data(connect=connect)
     #plt.axis('equal')
     self.ax.set_aspect('equal', 'datalim')
     if self.center is not None:
         cx,cy=self.center.real,self.center.imag; r=self.radius
         self.ax.axis([cx-w*r,cx+w*r,cy-w*r,cy+w*r])
     else:
         xmx=amax(self.x); ymn,ymx=amin(self.y),amax(self.y)
         cx=0.5*xmx; cy=0.5*(ymn+ymx); r=0.5*(ymx-ymn)
         self.ax.axis([cx-w*r,cx+w*r,cy-w*r,cy+w*r])
     if self.ZorY == 'Z':
         plt.xlabel(r'resistance $R$ in Ohm'); plt.ylabel(r'reactance $X$ in Ohm')
     if self.ZorY == 'Y':
         plt.xlabel(r'conductance $G$ in Siemens'); plt.ylabel(r'susceptance $B$ in Siemens')

def plotter(mode,Bc,Tc,Q):
    col = ['#000080','#0000FF','#4169E1','#6495ED','#00BFFF','#B0E0E6']
    plt.figure()
    ax = plt.subplot(111)
    for p in range(Bc.shape[1]):
        plt.plot(Tc[:,p],Bc[:,p],'-',color=str(col[p]))
    plt.xlabel('Tc [TW]')
    plt.ylabel('Bc normalised to total EU load')
    plt.title(str(mode)+' flow')
    
    # Shrink current axis by 25% to make room for legend
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width * 0.75, box.height])

    plt.legend(\
        ([str(Q[i]*100) for i in range(len(Q))]),\
        loc='center left', bbox_to_anchor=(1, 0.5),title='Quantiles')
    
    plt.savefig('figures/bctc_'+str(mode)+'.eps')

def _plot_base(dep, val, deplim_small, xlim_small, xlabel):
    plt.subplot(1,2,1)
    plt.plot(val, dep)
    plt.gca().invert_yaxis()
    plt.grid('on')
    plt.ylabel('depth/km')
    plt.xlabel(xlabel)
    locs, labels = plt.xticks()
    plt.setp(labels, rotation=-45)

    plt.subplot(1,2,2)
    plt.plot(val, dep)
    plt.gca().invert_yaxis()
    plt.grid('on')
    plt.ylim(deplim_small)
    plt.xlim(xlim_small)
    plt.xlabel(xlabel)
    locs, labels = plt.xticks()
    plt.setp(labels, rotation=-45)

def freqz(sosmat, nsamples=44100, sample_rate=44100, plot=True):
    """Plots Frequency response of sosmat."""
    from pylab import np, plt, fft, fftfreq
    x = np.zeros(nsamples)
    x[nsamples/2] = 0.999
    y, states = sosfilter_double_c(x, sosmat)
    Y = fft(y)
    f = fftfreq(len(x), 1.0/sample_rate)
    if plot:
        plt.grid(True)
        plt.axis([0, sample_rate / 2, -100, 5])
        L = 20*np.log10(np.abs(Y[:len(x)/2]) + 1e-17)
        plt.semilogx(f[:len(x)/2], L, lw=0.5)
        plt.hold(True)
        plt.title('freqz sos filter')
        plt.xlabel('Frequency / Hz')
        plt.ylabel('Damping /dB(FS)')
        plt.xlim((10, sample_rate/2))
        plt.hold(False)
    return x, y, f, Y