added optimization

Fa2h.py

@@ -0,0 +1,52 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Feb 04 17:59:00 2016
+from Matlab
+@author: Tobias Heinl, Gerald Schuller
+
+Function extracts analysis baseband impulse response (reverse window function) 
+from the folding matrix Fa of a cosine 
+modulated filter bank, with modulation function for the analysis IR : 
+h_k(n)=h(n)*cos(pi/N*(k+0.5)*(L-1-n+0.5-N/2));
+
+"""
+
+
+def Fa2h(Fa):
+   """ Function extracts analysis baseband impulse response (reverse window function) 
+   from the folding matrix Fa of a cosine modulated filter bank.
+   """
+   import numpy as np
+   from polmatmult import polmatmult
+   [N,y,blocks] = np.shape(Fa) 
+   h0 = np.zeros(blocks * N)
+   #First column of DCT-4:
+   T = np.zeros((N,1,1)) 
+   T[:,0,0]=np.cos(np.pi/N*(0.5)*(np.arange(N)+0.5))
+   #Compute first column of Polyphase matrix Pa(z):
+   Pa =  polmatmult(Fa,T)
+   #Extract impulse response h0(n):
+   for m in range(blocks): 
+      h0[m*N+np.arange(N)] = np.flipud(Pa[:,0,m]) 
+   #Baseband prototype h(n), divide by modulation func.:
+   #h = h0 / np.cos(np.pi/N*0.5*(np.arange(blocks*N)+ 0.5+(N/2)))
+   h = -h0 / np.cos(np.pi/N*0.5*(np.arange(blocks*N-1,-1,-1)+ 0.5-(N/2)))
+   return h;  
+
+
+#Testing:
+if __name__ == '__main__':
+  import numpy as np
+  from symFmatrix import *
+  from Dmatrix import *
+  from polmatmult import *
+  f=np.arange(1,7)
+  print("f=", f)
+  Fa=symFmatrix(f)
+  N=4
+  D=Dmatrix(N)
+  Faz=polmatmult(Fa,D)
+  print("Faz=", Faz[:,:,0],"\n", Faz[:,:,1])
+  h=Fa2h(Faz)
+  print("h=", h)
+

MDCTfb.py

@@ -57,9 +57,9 @@ def MDCTsynfb(y,fb):
    N=4
    D=Dmatrix(N)
    Dinv=Dinvmatrix(N)
-   #Filter bank coefficients, 1.5*N of sine window:
+   #Filter bank coefficients for sine window:
    #fb=np.sin(np.pi/(2*N)*(np.arange(int(1.5*N))+0.5))
-   fb=np.loadtxt("MDCTcoeff.txt")
+   fb=np.loadtxt("MDCTcoeff.txt") #Coeff. from optimization
    print("fb=", fb)
    #input test signal, ramp:
    x=np.arange(64)

PythonPsychoacoustics/psyac_quantization.py

@@ -23,7 +23,7 @@ def MDCT_psayac_quant_enc(x,fs,fb,N, nfilts=64,quality=100):
    #y: the MDCT subbands,
    #mTbarkquant: The quantied masking threshold in the Bark domain
    #usage: yq, y, mTbarkquant =MDCT_psayac_quant_enc(x,fs,fb,N)
-
+   print("quality=", quality)
    maxfreq=fs/2
    alpha=0.8  #Exponent for non-linear superposition of spreading functions
    nfft=2*N  #number of fft subbands
@@ -51,7 +51,7 @@ def MDCT_psayac_quant_enc(x,fs,fb,N, nfilts=64,quality=100):
    for m in range(M): #M: number of blocks
      #Observe: MDCT instead of STFT as input can be used by replacing ys by y:
      mXbark=mapping2bark(np.abs(ys[:,m]),W,nfft)
-     mTbark=maskingThresholdBark(mXbark,spreadingfuncmatrix,alpha,fs,nfilts)/(quality/100.0)
+     mTbark=maskingThresholdBark(mXbark,spreadingfuncmatrix,alpha,fs,nfilts)/(quality/100)
      #Logarithmic quantization of the "scalefactors":
      mTbarkquant[:,m]=np.round(np.log2(mTbark)*4) #quantized scalefactors
      mTbarkquant[:,m]=np.clip(mTbarkquant[:,m],0,None) #dequantized is at least 1 
@@ -127,7 +127,7 @@ def MDCTsyn_dequant_dec(yq, mTbarkquant, fs, fb, N, nfilts=64):
    fb=np.sin(np.pi/(2*N)*(np.arange(int(1.5*N))+0.5))
    print("Encoder part:")
    #MDCT and quantization:
-   yq, y, mTbarkquant = MDCT_psayac_quant_enc(x,fs,fb,N, nfilts,quality=60.0)
+   yq, y, mTbarkquant = MDCT_psayac_quant_enc(x,fs,fb,N, nfilts,quality=60)
 
    print("Decoder part:")
    xrek, mT, ydeq = MDCTsyn_dequant_dec(yq, mTbarkquant, fs, fb, N, nfilts)
@@ -172,7 +172,7 @@ def MDCTsyn_dequant_dec(yq, mTbarkquant, fs, fb, N, nfilts=64):
    plt.ylabel("Normalized Frequency (pi is Nyquist freq.)")
    plt.show()
 
-   plt.plot(xrek)
+   plt.plot(x)
    plt.title("The reconstructed audio signal")
    plt.show()
 

PythonPsychoacoustics/psyacmodel.py

@@ -24,6 +24,7 @@ def f_SP_dB(maxfreq,nfilts):
 def spreadingfunctionmat(spreadingfunctionBarkdB,alpha,nfilts):
    #Turns the spreading prototype function into a matrix of shifted versions.
    #Convert from dB to "voltage" and include alpha exponent
+   #nfilts: Number of subbands in the Bark domain, for instance 64  
    spreadingfunctionBarkVoltage=10.0**(spreadingfunctionBarkdB/20.0*alpha)
    #Spreading functions for all bark scale bands in a matrix:
    spreadingfuncmatrix=np.zeros((nfilts,nfilts))
@@ -33,17 +34,16 @@ def spreadingfunctionmat(spreadingfunctionBarkdB,alpha,nfilts):
 
 
 
-
 def maskingThresholdBark(mXbark,spreadingfuncmatrix,alpha,fs,nfilts): 
   #Computes the masking threshold on the Bark scale with non-linear superposition
   #usage: mTbark=maskingThresholdBark(mXbark,spreadingfuncmatrix,alpha)
   #Arg: mXbark: magnitude of FFT spectrum, on the Bark scale
   #spreadingfuncmatrix: spreading function matrix from function spreadingfunctionmat
   #alpha: exponent for non-linear superposition (eg. 0.6), 
   #fs: sampling freq., nfilts: number of Bark subbands
-  #return: masking threshold as "voltage" on Bark scale 
-  #Returns:
-  #mTbark: the resulting Masking Threshold on the Bark scale 
+  #nfilts: Number of subbands in the Bark domain, for instance 64  
+  #Returns: mTbark: the resulting Masking Threshold on the Bark scale 
+
   #Compute the non-linear superposition:
   mTbark=np.dot(mXbark**alpha, spreadingfuncmatrix**alpha)
   #apply the inverse exponent to the result:

README.md

@@ -55,6 +55,21 @@ For the binary file, the library "pickle" is needed, it can be installed with
 
 sudo pip3 install pickle 
 
+** Filter Bank Optimization
 
-Gerald Schuller, gerald.schuller@tu-ilmenau.de, June 2018.
+The folder also contains a few programs which show how to optimize different types of filter banks, with regard to their filter characteristics.
+
+The following runs an optimization for an MDCT, in this example for N=4 subbands (and filter length 2N=8),
+
+python3 optimfuncMDCT.py
+
+This runs an optimization for a Low Delay Filter Bank, also for N=4, but filter length 3N=12, and system delay of 7 (including blocking delay of N-1=3, which doesn't show up in file based examples),
+
+python3 optimfuncLDFB.py
+
+The last runs an optimization for a PQMF filter bank, for N=4 subbands and filter length 8N=32, and system delay of 31 (including the blocking delay of N-1=3),
+
+python3 optimfuncQMF.py
+
+Gerald Schuller, gerald.schuller@tu-ilmenau.de, July 2018.
 

optimfuncLDFB.py

@@ -0,0 +1,77 @@
+# -*- coding: utf-8 -*-
+"""
+Program to optimize cosine modulated low delay filter banks.
+
+Gerald Schuller, July 2018
+"""
+
+import numpy as np
+import scipy as sp
+import scipy.signal as sig
+from symFmatrix import symFmatrix
+from polmatmult import polmatmult
+from Dmatrix import Dmatrix
+from Fa2h import Fa2h
+from Gmatrix import Gmatrix
+
+def optimfuncLDFB(x,N):
+    '''function for optimizing an MDCT type filter bank.
+    x: unknown matrix coefficients, N: Number of subbands.
+    '''
+    #Analysis folding matrix:
+    Fa = symFmatrix(x[0:int(1.5*N)])
+    Faz = polmatmult(Fa,Dmatrix(N))
+    Faz = polmatmult(Faz,Gmatrix(x[int(1.5*N):(2*N)]))
+    #baseband prototype function h:
+    h = Fa2h(Faz)
+    #'Fa2h is returning 2D array. Squeeze to make it 1D'
+    #h= np.squeeze(h)
+
+    #Frequency response of the the baseband prototype function at 1024 frequency sampling points
+    #between 0 and Nyquist:
+    w,H = sig.freqz(h,1,1024)
+    #desired frequency response
+    #Limit of desired pass band (passband is between -pb and +pb, hence ../2)
+    pb = int(1024/N/2.0)
+    #Ideal desired frequency response:
+    Hdes = np.concatenate((np.ones(pb),np.zeros(1024-pb)))
+    #transition band width to allow filter transition from pass band to stop band:
+    tb = int(np.round(1.0*pb))
+    #Weights for differently weighting errors in pass band, transition band and stop band:
+    weights = np.concatenate((1.0*np.ones(pb), np.zeros(tb),1000*np.ones(1024-pb-tb)))
+    #Resulting total error number as the sum of all weighted errors:
+    err = np.sum(np.abs(H-Hdes)*weights)
+    return err
+
+if __name__ == '__main__':  #run the optimization
+  import scipy as sp
+  import scipy.signal
+  import matplotlib.pyplot as plt
+
+  N=4  #number of subbands
+  s=2*N
+  bounds=[(-14,14)]*s
+  xmin = sp.optimize.differential_evolution(optimfuncLDFB, bounds, args=(N,), disp=True)
+  print("error after optim.=",xmin.fun)
+  print("optimized coefficients=",xmin.x)
+  np.savetxt("LDFBcoeff.txt", xmin.x)
+  x=xmin.x;
+  #Baseband Impulse Response:
+  Fa = symFmatrix(x[0:int(1.5*N)])
+  #print("Fa=", Fa[:,:,0])
+  Faz = polmatmult(Fa,Dmatrix(N))
+  Faz = polmatmult(Faz,Gmatrix(x[int(1.5*N):(2*N)]))
+  h = Fa2h(Faz)
+  plt.plot(h)
+  plt.xlabel('Sample')
+  plt.ylabel('Value')
+  plt.title('Baseband Impulse Response of the Low Delay Filter Bank')
+  #Magnitude Response:
+  w,H=sp.signal.freqz(h) 
+  plt.figure()
+  plt.plot(w,20*np.log10(abs(H)))
+  plt.axis([0, 3.14, -60,20])
+  plt.xlabel('Normalized Frequency')
+  plt.ylabel('Magnitude (dB)')
+  plt.title('Mag. Frequency Response of the Low Delay Filter Bank')
+  plt.show()

optimfuncMDCT.py

@@ -0,0 +1,69 @@
+# -*- coding: utf-8 -*-
+"""
+Program for optimizing an MDCT type filter bank with N subbands. 
+Gerald Schuller, July 2018
+"""
+
+
+def optimfuncMDCT(x, N):
+    """Computes the error function for the filter bank optimization
+    for coefficients x, a 1-d array, N: Number of subbands""" 
+    import numpy as np
+    import scipy.signal as sig
+    from polmatmult import polmatmult 
+    from Dmatrix import Dmatrix
+    from symFmatrix import symFmatrix
+    from Fa2h import Fa2h
+
+    #x = np.transpose(x) 
+    Fa = symFmatrix(x)
+    D = Dmatrix(N)
+    Faz = polmatmult(Fa,D)
+    h = Fa2h(Faz)
+    h = np.hstack(h)
+    w, H = sig.freqz(h,1,1024)
+    pb = int(1024/N/2)
+    Hdes = np.concatenate((np.ones((pb,1)) , np.zeros(((1024-pb, 1)))), axis = 0)
+    tb = np.round(pb)
+    weights = np.concatenate((np.ones((pb,1)) , np.zeros((tb, 1)), 1000*np.ones((1024-pb-tb,1))), axis = 0) 
+    err = np.sum(np.abs(H-Hdes)*weights)
+    return err  
+
+if __name__ == '__main__':  #run the optimization
+  import numpy as np
+  import scipy as sp
+  import scipy.optimize
+  import scipy.signal
+  import matplotlib.pyplot as plt
+  from symFmatrix import symFmatrix
+  from polmatmult import polmatmult 
+  from Dmatrix import Dmatrix
+  from Fa2h import Fa2h
+
+  N=4
+  #Start optimization with some starting point:
+  x0 = -np.random.rand(int(1.5*N))
+  print("starting error=", optimfuncMDCT(x0, N)) #test optim. function
+  xmin = sp.optimize.minimize(optimfuncMDCT, x0, args=(N,), options={'disp':True})
+  print("optimized coefficients=", xmin.x)
+  np.savetxt("MDCTcoeff.txt", xmin.x)
+  print("error after optim.=", optimfuncMDCT(xmin.x, N))
+  #Baseband Impulse Response:
+  Fa = symFmatrix(xmin.x) 
+  Faz = polmatmult(Fa, Dmatrix(N))
+  h = Fa2h(Faz)
+  print("h=", h)
+  plt.plot(h)
+  plt.xlabel('Sample')
+  plt.ylabel('Value')
+  plt.title('Baseband Impulse Response of our Optimized MDCT Filter Bank')
+  plt.figure()
+  #Magnitude Response:
+  w,H=sp.signal.freqz(h) 
+  plt.plot(w,20*np.log10(abs(H)))
+  plt.axis([0, 3.14, -60,20])
+  plt.xlabel('Normalized Frequency')
+  plt.ylabel('Magnitude (dB)')
+  plt.title('Mag. Frequency Response of the MDCT Filter Bank') 
+  plt.show()   
+

optimfuncQMF.py

@@ -0,0 +1,72 @@
+# -*- coding: utf-8 -*-
+"""
+Python programs to optimize an PQMF filter bank with N subbands
+Gerald Schuller, July 2018
+"""
+
+
+from __future__ import print_function
+def optimfuncQMF(x,N):
+    """Optimization function for a PQMF Filterbank
+    x: coefficients to optimize (first half of prototype h), N: Number of subbands
+    err: resulting total error
+    """
+    import numpy as np
+    import scipy as sp
+    import scipy.signal as sig 
+
+    h = np.append(x,np.flipud(x));
+    #H = sp.freqz(h,1,512,whole=True)
+    f,H_im = sig.freqz(h)
+    H=np.abs(H_im) #only keeping the real part
+    posfreq = np.square(H[0:int(512/N)]);
+    #Negative frequencies are symmetric around 0:
+    negfreq = np.flipud(np.square(H[0:int(512/N)]))
+    #Sum of magnitude squared frequency responses should be close to unity (or N)
+    unitycond = np.sum(np.abs(posfreq+negfreq - 2*(N*N)*np.ones(int(512/N))))/512;
+    #plt.plot(posfreq+negfreq);
+    #High attenuation after the next subband:
+    att = np.sum(np.abs(H[int(1.5*512/N):]))/512;
+    #Total (weighted) error:
+    err = unitycond + 100*att;
+    return err
+
+if __name__ == '__main__':  #run the optimization
+   from scipy.optimize import minimize
+   import scipy as sp
+   import matplotlib.pyplot as plt
+
+   N=4 #Number of subbands
+   #Start optimization with "good" starting point:
+   x0 = 16*sp.ones(4*N)
+   print("starting error=", optimfuncQMF(x0,N)) #test optim. function
+   xmin = minimize(optimfuncQMF,x0, args=(N,), method='SLSQP')
+   print("error after optim.=",xmin.fun)
+   print("optimized coefficients=",xmin.x)
+   #Store the found coefficients in a text file: 
+   sp.savetxt("QMFcoeff.txt", xmin.x)
+   #we compute the resulting baseband prototype function:
+   h = sp.concatenate((xmin.x,sp.flipud(xmin.x)))
+   plt.plot(h)
+   plt.xlabel('Sample')
+   plt.ylabel('Value')
+   plt.title('Baseband Impulse Response of the Optimized PQMF Filter Bank')
+   #plt.xlim((0,31))
+   plt.show()
+   #The corresponding frequency response:
+   w,H=sp.signal.freqz(h) 
+   plt.plot(w,20*sp.log10(abs(H)))
+   plt.axis([0, 3.14, -100,20])
+   plt.xlabel('Normalized Frequency')
+   plt.ylabel('Magnitude (dB)')
+   plt.title('Mag. Frequency Response of the PQMF Filter Bank')
+   plt.show()
+   #Checking the "unity condition":
+   posfreq = sp.square(abs(H[0:int(512/N)]));
+   negfreq = sp.flipud(sp.square(abs(H[0:int(512/N)])))
+   plt.plot(posfreq+negfreq)
+   plt.xlabel('Frequency (512 is Nyquist)')
+   plt.ylabel('Magnitude')
+   plt.title('Unity Condition, Sum of Squared Magnitude of 2 Neigh. Subbands')
+   plt.show() 
+