#!/usr/bin/python
import getopt, sys

from cvxopt import matrix
from cvxopt.solvers import options
from cvxopt.modeling import op, variable, max

from math import cos, log10, pi

import pygtk
pygtk.require('2.0')
import gtk

import matplotlib
matplotlib.use('GTKAgg')  # or 'GTK'
from matplotlib.backends.backend_gtk import FigureCanvasGTK as FigureCanvas
import pylab 

def frange(a,b,N):    
    return [ a+k*float((b-a))/N  for k in xrange(N) ]

def design_lowpass(N, d1, wc, ws, solver=None, Q=50):

    h = variable(N+1)
    d1 = 10**(d1/20.0)     # convert from dB
    
    n1 = int(round(N*Q*wc/pi));
    w1 = matrix(frange(0,wc,n1))
    G1 = matrix([cos(wi*j) for j in xrange(N+1) for wi in w1], (n1,N+1))

    n2 = int(round(N*Q*(pi-ws)/pi));
    w2 = matrix(frange(ws,pi,n2))
    G2 = matrix([cos(wi*j) for j in xrange(N+1) for wi in w2], (n2,N+1))
    
    options['show_progress'] = 0    
    options['LPX_K_MSGLEV'] = 0
    options['MSK_IPAR_LOG']= 0
    op(max(abs(G2*h)), [G1*h <= d1, G1*h >= 1.0/d1]).solve(solver=solver)

    return (h.value, max(abs(G2*h.value)))

def make_plot(h, d2, co, sb, pr, N, output=None):
    w = matrix(frange(0,pi,N*50));        
    C = w*matrix(range(N+1), (1,N+1), 'd');
    for i in xrange(len(C)): C[i] = cos(C[i])
    
    fig = pylab.figure()
    ax = fig.add_subplot(111)
    ylim = [round((20*log10(d2)-40)/10)*10,10]
    ax.plot(list(w/pi),[20*log10(abs(x)) for x in C*h])
    ax.plot(2*[co],[-pr,ylim[0]],'g--')
    ax.plot(2*[co],[10,pr],'g--')
    ax.plot(2*[sb],[10,20*log10(d2)],'g--')
    ax.plot([sb, 1],2*[20*log10(d2)],'g--')
    ax.plot([0, co],2*[-pr],'g--')
    ax.plot([0, co],2*[pr],'g--')

    pylab.setp(ax,ylim=ylim)
    pylab.setp(ax,xlabel="Normalized frequency")
    pylab.setp(ax,ylabel="Attenuation [dB]")
    ax.grid()

def usage():
    print """
Usage:
filterdemo_cli --cutoff=CO --stopband=SB --ripple=RP --order=N [options]

Arguments:
 CO: normalized cutoff frequency
 SB: normalized stopband frequency, 0.1 <= co < sb-0.1 <= 0.5,
 RP: maximum passband ripple in dB, 0.01 <= rp <= 3,
 N : filterorder, 5 <= or <= 50.

Options:
--solver = SOLVER      One of default, mosek, glpk
--output = FILENAME    Output filename.
"""
    sys.exit(2)
    
def main():
    try:
        opts, args = getopt.getopt(sys.argv[1:], "",
             ['cutoff=', 'stopband=', 'ripple=', 'order=',
              'solver=', 'output='])
    except getopt.GetoptError:
        usage()

    if opts==[]: usage()
        
    co, sb, pr, N, output = [None]*5
    solver = "default"
    
    try:
        for o, a in opts:
            if o == "--cutoff":   co = float(a);
            if o == "--stopband": sb = float(a);
            if o == "--ripple":   pr = float(a);
            if o == "--order":    N  = int(a);
            if o == "--solver":   solver = a;
            if o == "--output":   output = a;            
    except: usage()

    if None in [co, sb, pr, N]: usage()
    
    if not (0.1 <= co < sb-0.01+1E-8 <= 0.5):
        print "invalid cutoff and stopband frequencies"
        usage()

    if not (0.01 <= pr <= 3):
        print "invalid value of passband ripple"
        usage()

    if not (5 <= N <= 50):
        print "invalid filterorder"
        usage()
    
    if not solver in ['default','mosek','glpk']:
        print "invalid solver"
        usage()


    try:
        [h, d2] = design_lowpass(N, pr, co*pi, sb*pi, solver)
    except:
        print "Please tighten filter specifications."
        sys.exit(2)
        
        
    make_plot(h, d2, co, sb, pr, N, output);

    if (output != None): savefig(output)

    pylab.show()
        
if __name__ == "__main__":
    main()