import numpy as np
import matplotlib.pyplot as plt
import copy
from   scipy import signal as sig

# Generate anti aliasing filter
def genAAFilter(Att, M, factor):
    def Ino(x):
        # Zero order Bessel function
        d = 0
        ds = 1
        s = 1

        while(ds > s*1e-6):
            d += 2
            ds *= x*x/(d*d)
            s += ds

        return s;

    Fa = 0;
    Fb = 1 / (2*factor);

    Np = int((M-1)/2)
    A = np.zeros(M)
    A[0] = 2*(Fb-Fa);
    for j in range(1,Np+1):
        A[j] = (np.sin(2*j*np.pi*Fb)-np.sin(2*j*np.pi*Fa))/(j*np.pi);

    if (Att<21):
      Alpha = 0;
    elif (Att>50):
      Alpha = 0.1102*(Att-8.7);
    else:
      Alpha = 0.5842*Math.pow((Att-21), 0.4)+0.07886*(Att-21);

    # Window the ideal response with the Kaiser-Bessel window
    Inoalpha = Ino(Alpha);
    H = np.zeros(M)
    for j in range(Np+1):
      H[Np+j] = A[j]*Ino(Alpha*np.sqrt(1-(j*j/(Np*Np))))/Inoalpha;

    for j in range(Np+1):
      H[j] = H[M-1-j];

    return H

def genFilterCode(H, idx, skipFactor):
    string = '';
    t = 0;

    # Find unique weights (there could be duplicates for symmetry reasons)
    # For the windowed sinc, (N-1)/2 duplicates actually.
    weightRefs = [];
    for i in range(len(H)):
        nidx = i;
        for j in range(i):
            if H[i] == H[j]:
                nidx = j;
                break;
        weightRefs.append(nidx)

    wHist = ['']*len(H);
    # Assign history to weight list
    for i in range(idx, len(H), skipFactor):
        # Loop over the history bit that is active for this polyphase filter
        # t holds the current history index, while i holds the current weight index

        # Look up the unique-ified weight for this
        wIdx = weightRefs[i]

        # Assign the history element to this list
        wHist[wIdx] = wHist[wIdx] + ' + h%d'%t

        t = t + 1;

    # Start printing code
    for i in range(len(H)):
        if np.abs(H[i]) > 0:#1e-7:
            if len(wHist[i]) > 0:
                if H[i] > 0:
                    string = "%s + %s * (%s)" % (string, '{:.28f}'.format(H[i]), wHist[i][3:])
                    #string = "%s + w%d * (%s)" % (string, i, wHist[i][3:])
                else:
                    string = "%s - %s * (%s)" % (string, '{:.28f}'.format(abs(H[i])), wHist[i][3:])
                    #string = "%s - w%d * (%s)" % (string, i, wHist[i][3:])

    return string, t+1

def generateHeader(name, M, extra = ''):
    string = 'function %s\n' % name
    string = string + '  instance('
    for i in range(M):
        string = '%sh%d, ' % (string,i)
    string = string[0:-2];
    string = string + extra + ')\n  global()\n  local()\n';

    return string

def generateDownSampleFilter(name, M, factor, att):
    code = generateHeader(name, M)
    H = genAAFilter(att, M, factor);
    downSampleCode, maxHist = genFilterCode(H, 0, 1)
    code = code + '  (\n    ' + downSampleCode + '\n  );\n'

    return code, H, maxHist

def generateUpSampleFilter(name, M, factor, att):
    H = genAAFilter(att, M, factor);
    code = generateHeader(name, M, ", idx ")
    code = '%s  (\n    idx += 1;\n    ( idx == 0 ) ? (\n' % (code)

    for j in range(factor):
        upSampleCode, maxHist = genFilterCode(H, j, factor)
        if j < (factor-1):
            ifs =  ' : ( idx == %d ) ? (' % (j + 1)
        else:
            ifs = '';
        code = '%s      %s\n    )%s\n' % (code, upSampleCode, ifs)

    code = code + '  );';

    return code, H, maxHist

def generateInitializer(name, M):
    code = generateHeader(name, M)
    code = '%s  (\n     '%code

    for i in range(M):
        code = '%sh%d = ' % (code,i)

    code = '%s0\n  );' % code;

    return code;

def generateHistoryUpdater(name, M):
    code = generateHeader(name, M, ', idx')
    code = '%s  (\n'%code

    for i in range(1,M):
        code = '%s     h%d = h%d;\n' % (code,M-i,M-i-1)

    code = '%s\n     h0 = x;\n     idx = -1;\n  );' % code;

    return code;

# Attenuation in dB
Att = 96

fid  = open('saike_upsamplers.jsfx-inc', 'w')

# Filter lengths for different oversampling ratios
lens = [65, 97, 129, 161, 193, 225, 257];

# Print out delays, make sure that these aren't fractional.
print( (.5*(np.array(lens)-1)+.5*(np.array(lens)-1))/np.array([2,3,4,5,6,7,8]) )

fid.write('@init\n')
j = 0;
upList = [];

# Compute oversampling kernels from x2 to x8
for factor in np.arange(2,9):
    downCode, H, maxHistDown = generateDownSampleFilter( 'downSample%d()'%factor, lens[j], factor, 96 )
    upCode, H, maxHistUp = generateUpSampleFilter( 'upSample%d()'%factor, lens[j], factor, 96 )

    upList.append(maxHistUp);

    fid.write(generateInitializer('init_upSampler%d()'%factor, lens[j]));
    fid.write('\n\n')

    fid.write(generateHistoryUpdater('update_upHist%d(x)'%factor, maxHistUp));
    fid.write('\n\n')

    fid.write(generateHistoryUpdater('update_downHist%d(x)'%factor, maxHistDown));
    fid.write('\n\n')

    fid.write(downCode)
    fid.write('\n\n')
    fid.write(upCode)
    fid.write('\n\n')
    j = j + 1;

fid.write('function upSample(N)\n')
fid.write('  instance()\n')
fid.write('  global()\n')
fid.write('  local()\n')
fid.write('(')
fid.write(' (N==2) ? this.upSample2()\n')
for factor in np.arange(3,9):
    fid.write(' : (N==%d) ? this.upSample%d()\n'%(factor,factor))
fid.write(');\n\n')
fid.write('function downSample(N)\n')
fid.write('  instance()\n')
fid.write('  global()\n')
fid.write('  local()\n')
fid.write('(\n')
fid.write(' (N==2) ? this.downSample2()\n')
for factor in np.arange(3,9):
    fid.write(' : (N==%d) ? this.downSample%d()\n'%(factor,factor))
fid.write(');\n\n')
fid.write('function updateDownHist(N, x)\n')
fid.write('  instance()\n')
fid.write('  global()\n')
fid.write('  local()\n')
fid.write('(\n')
fid.write(' (N==2) ? this.update_downHist2(x)\n')
for factor in np.arange(3,9):
    fid.write(' : (N==%d) ? this.update_downHist%d(x)\n'%(factor,factor))
fid.write(');\n\n')
fid.write('function updateUpHist(N, x)\n')
fid.write('  instance()\n')
fid.write('  global()\n')
fid.write('  local()\n')
fid.write('(\n')
fid.write(' (N==2) ? this.update_upHist2(x)\n')
for factor in np.arange(3,9):
    fid.write(' : (N==%d) ? this.update_upHist%d(x)\n'%(factor,factor))
fid.write(');\n\n')

fid.write('function getFIRdelay(N)\n')
fid.write('  instance()\n')
fid.write('  global()\n')
fid.write('  local()\n')
fid.write('(\n')
fid.write('    (N==1) ? 0\n');
for j in range(len(upList)):
    fid.write(' : (N==%d) ? %d\n'%(j+2,(upList[j]-1)))
fid.write(');\n\n')


fid.close()