TTS/vocoder/notebooks/Untitled.ipynb

679 lines
226 KiB
Plaintext
Raw Normal View History

{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"#function example with several unknowns (variables) for optimization\n",
"#Gerald Schuller, Nov. 2016\n",
"import numpy as np\n",
"\n",
"def functionexamp(x):\n",
" #x: array with 2 variables\n",
" \n",
" y=np.sin(x[0])+np.cos(x[1])\n",
" return y"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" fun: -1.9999999999888387\n",
" jac: array([4.7236681e-06, 0.0000000e+00])\n",
" message: 'Optimization terminated successfully.'\n",
" nfev: 12\n",
" nit: 2\n",
" njev: 3\n",
" status: 0\n",
" success: True\n",
" x: array([-1.5707916 , -3.14159265])\n"
]
}
],
"source": [
"#Optimization example, see also:\n",
"#https://docs.scipy.org/doc/scipy-0.18.1/reference/optimize.html\n",
"#Gerald Schuller, Nov. 2016\n",
"#run it with \"python optimizationExample.py\" in a termina shell\n",
"#or type \"ipython\" in a termina shell and copy lines below:\n",
"\n",
"import numpy as np\n",
"import scipy.optimize as optimize\n",
"\n",
"#Example for 2 unknowns, args: function-name, starting point, method:\n",
"xmin = optimize.minimize(functionexamp, [-1.0, -3.0], method='CG')\n",
"print(xmin)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"function [p,passedge] = opt_filter(filtorder,N)\n",
"\n",
"% opt_filter Create Lowpass Prototype Filter for the Pseudo-QMF \n",
"% Filter Bank with N Subbands\n",
"%\n",
"% Adapted from the paper by C. D. Creusere and S. K. Mitra, titled \n",
"% \"A simple method for designing high-quality prototype filters for \n",
"% M-band pseudo-QMF banks,\" IEEE Trans. Signal Processing,vol. 43, \n",
"% pp. 1005-1007, Apr. 1995 and the book by S. K. Mitra titled \"\n",
"% Digital Signal Processing: A Computer-Based Approach, McGraw-Hill, 2001\n",
"%\n",
"% Arguments:\n",
"% filtorder Filter order (i.e., filter length - 1)\n",
"% N Number of subbands\n",
"\n",
"stopedge = 1/N; % Stopband edge fixed at (1/N)pi\n",
"passedge = 1/(4*N); % Start value for passband edge\n",
"tol = 0.000001; % Tolerance\n",
"step = 0.1*passedge; % Step size for searching the passband edge\n",
"way = -1; % Search direction, increase or reduce the passband edge\n",
"tcost = 0; % Current error calculated with the cost function\n",
"pcost = 10; % Previous error calculated with the cost function\n",
"flag = 0; % Set to 1 to stop the search\n",
"\n",
"while flag == 0\n",
" \n",
"% Design the lowpass filter using Parks-McClellan algorithm\n",
" \n",
" p = remez(filtorder,[0,passedge,stopedge,1],[1,1,0,0],[5,1]);\n",
" \n",
"% Calculates the cost function according to Eq. (2.36)\n",
"\n",
" P = fft(p,4096);\n",
" OptRange = floor(2048/N); % 0 to pi/N\n",
" phi = zeros(OptRange,1); % Initialize to zeros\n",
"\n",
"% Compute the flatness in the range from 0 to pi/N\n",
"\n",
"\tfor k = 1:OptRange\n",
" phi(k) = abs(P(OptRange-k+2))^2 + abs(P(k))^2;\n",
"\tend\n",
"\ttcost = max(abs(phi - ones(max(size(phi)),1)));\n",
" \t\n",
"\tif tcost > pcost % If search in wrong direction\n",
"\t\tstep = step/2; % Reduce step size by half \n",
"\t\tway = -way; % Change the search direction \n",
"\tend\n",
"\t\n",
"\tif abs(pcost - tcost) < tol % If improvement is below tol \n",
"\t\tflag = 1; % Stop the search \n",
"\tend\n",
"\t\n",
"\tpcost = tcost;\n",
"\tpassedge = passedge + way*step; % Adjust the passband edge\n",
" \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"sig.remez"
]
},
{
"cell_type": "code",
"execution_count": 101,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"text/plain": [
"0.0125"
]
},
"execution_count": 101,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"1 / 4. / 20.0"
]
},
{
"cell_type": "code",
"execution_count": 90,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"ename": "ValueError",
"evalue": "Band edges should be less than 1/2 the sampling frequency",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-90-fb4445f5709b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mremez\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m64\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m16.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m4.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m~/miniconda3/lib/python3.7/site-packages/scipy/signal/fir_filter_design.py\u001b[0m in \u001b[0;36mremez\u001b[0;34m(numtaps, bands, desired, weight, Hz, type, maxiter, grid_density, fs)\u001b[0m\n\u001b[1;32m 854\u001b[0m \u001b[0mbands\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbands\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 855\u001b[0m return sigtools._remez(numtaps, bands, desired, weight, tnum, fs,\n\u001b[0;32m--> 856\u001b[0;31m maxiter, grid_density)\n\u001b[0m\u001b[1;32m 857\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 858\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mValueError\u001b[0m: Band edges should be less than 1/2 the sampling frequency"
]
}
],
"source": [
"p = sig.remez(65, [0, 1/16.0, 1/4.0, 1], [1, 0], [5, 1])\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"def create_pqmf_filter(filter_len=64, N=4):\n",
" stop_edge = 1 / N\n",
" pass_edge = 1 / (4 * N)\n",
" tol = 1e-8\n",
" cutoff = 0.1 * pass_edge\n",
" cost = 0\n",
" cost_prev = float('inf')\n",
" \n",
" p = sig.remez(filter_len, [0, pass_edge, stop_edge, 1], [1, 1, 0, 0], [5, 1])\n",
" \n",
" P = sig.freqz(p, workN=2048)\n",
" opt_range = 2048 // N\n",
" phi = np.zeros(opt_range)\n",
" \n",
" H = np.abs(P)\n",
" phi = H[opt_range + 2] \n",
" for i in range(opt_range):\n",
" phi[i] = abs(P(opt_range - i + 2)) ** 2 + abs(P[i]) ** 2"
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy as sp\n",
"import scipy.signal as sig\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"\n",
"\n",
"def optimfuncQMF(x):\n",
" \"\"\"Optimization function for a PQMF Filterbank\n",
" x: coefficients to optimize (first half of prototype h because of symmetry)\n",
" err: resulting total error\n",
" \"\"\"\n",
" K = ntaps * N \n",
" h = np.append(x, np.flipud(x))\n",
" cutoff = 0.15\n",
" \n",
"# breakpoint()\n",
" f, H_im = sig.freqz(h, worN=K)\n",
" H = np.abs(H_im) #only keeping the real part\n",
" \n",
" posfreq = np.square(H[0:K//N])\n",
" \n",
" #Negative frequencies are symmetric around 0:\n",
" negfreq = np.flipud(np.square(H[0:K//N]))\n",
" \n",
" #Sum of magnitude squared frequency responses should be closed to unity (or N)\n",
" unitycond = np.sum(np.abs(posfreq + negfreq - 2*(N*N)*np.ones(K//N)))/K\n",
" \n",
" #plt.plot(posfreq+negfreq)\n",
" \n",
" #High attenuation after the next subband:\n",
" att = np.sum(np.abs(H[int(cutoff*K//N):]))/K\n",
" \n",
" #Total (weighted) error:\n",
" err = unitycond + 100*att\n",
" return err"
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"text/plain": [
"(32,)"
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xmin.shape"
]
},
{
"cell_type": "code",
"execution_count": 86,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"8.684549400499243\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEWCAYAAABIVsEJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADt0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjByYzMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy9h23ruAAAgAElEQVR4nO3dd3xUVfr48c+TRgglhSQQEkLvHQMIoqCAYEVX17Wsoqvr6vbid7/6dVdd19113aJbv/74KnaxF0SFRQTpJbQQauikkUYSEtLz/P6YC44xIT035Xm/XvPK3Dvnzjz3zmSeOefce46oKsYYY0x9+bgdgDHGmLbJEogxxpgGsQRijDGmQSyBGGOMaRBLIMYYYxrEEogxxpgGsQRiWpyIrBKRe87z+LMi8uuWjKmtEJFYESkQEd8Gbn9URGY1dVytiYjsFpEZbsfREVgC6eCcL5Qi50spXUReFJGuLfj6d4rIWu91qnqfqv62GV6rn4ios68Fzr4/6PW4iMh/iUiSc0yOi8jvRSTAq8yLznPMq/LcTzvr7/Tarwqv1yoQkX9WE9MtIrK3yrrlNax7UFWPq2pXVa1oosPSpETkMREpq7Lfv2zG13tRRJ7wXqeqI1V1VXO9pvmSJRADcI2qdgXGAeOBh1yOp7mFOPt7C/CIiMx11v8duBe4A+gGXAHMAt6osv0BpwwAIuIH3AQcqlJug/Nlf/b2w2piWQ0ME5EIr+caC3Susm6KU7YteLPKfj/ldkCmeVgCMeeoajqwDE8iAUBELhSR9SKSKyI7vZsGnF/Zh0XktIgcEZHbnPWPicirXuXO/vL38349ERkOPAtMcX6p5jrrz/2qFJEZIpIsIr8QkQwRSRORu7yeo4eIfCQi+SKyRUSeqFqjOc/+bgB2A6NEZDDwfeA2Vd2gquWquhu4AbhKRKZ7bfoRME1EQp3luUACkF6X160SQwpwGLjEWTXBiemLKut8gC1Vj6XTHPhbEVnnvA//EZFwr+Nzu4gcE5FsEXnY+7VFpJOIPCMiqc7tGRHp5Dz2hYjc4Ny/yHnNq5zlmSKyoz77Wdtnog77Mc3rc3jC+ezdC9wG/NL5/HzklD3XTFfLPp73s2VqZwnEnCMiMXh+dR90lqOBj4EngDDgAeBdEYkQkS54frFfoardgKlAvb5UVHUvcB9f/lIPqaFoLyAYiAbuBv7l9eX9L6DQKTPfudVlX0VELgJGAtuBmUCyqm6uEuMJYCNwudfqYuBD4GZn+Q7g5bq8bg1W82WyuARYA6ytsm6jqpbVsP2twF1AJBCA531CREYA/wvcDvQGegAxXts9DFyI5wfDWGAS8CvnsS+AGc796Xw1yU13Hm9qNe1HX+BT4B9AhBPvDlVdALwGPOV8fq6p5jnPt49w/s+WqYUlEAPwgYicBk4AGcCjzvpvA5+o6ieqWqmqy4F44Ern8Uo8v947q2qa84u9OZQBj6tqmap+AhQAQ8XTkXwD8KiqnlHVPcBLdXi+LCAHeA54UFVXAOFAWg3l0/B8cXl7GbhDRELwfKF+UM12Fzq/mM/eLqzh+b1rGxfjSSBrqqw73xf2C6p6QFWLgLf4sgZ5I7BEVVeragnwazzv2Vm34TmuGaqaCfwGT7I5G9PZWtclwB+8lmtLIDdV2e/e5ylbl/24FfhMVRc5n4FsVa3rj5Xz7SPU8Nmq43N3eJZADMB1Ti1iBjAMz5cpQF/gm95fBsA0IEpVC4Fv4alBpInIxyIyrJniy1bVcq/lM0BXPF/qfngS31ne92sSrqqhqjpcVf/urMsComooH+U8fo6qrnVe/2E8X9JF1Wy3UVVDvG4ba3j+1cAY55fvhXhqZPuAKGfdNM7f/+HddHb22ICn1nHueDjvWbZX2d7AMa/lY846gA3AEBHpieeL/GWgj9OsNKmWeN6qst+p5ylbl/3ow9f7l+rqfPsINX+2TB1YAjHnqOoXwIvAn51VJ4BXqnwZdFHVJ53yy1R1Np4v2H3A/znbFQJBXk/d63wv24iQM4Fyvtos06eBz/U5ni/ISd4rRaQPni/1VdVs8yrwCxrXfIWqHgZS8XTgH1fVAuehDc66rnia0eorDa/jISJBeJqxzkrF8yPhrFhnHap6BtgK/ARIVNVSYD3wc+CQqn4lodZBfT4TVZ0ABtbwWG2fnxr30TSeJRBT1TPAbBEZi+cL8hoRmSMiviIS6HQ8xohITxGZ5/SFlOCp+p9tHtkBXCKeaxaCOf9ZXSeBGPE6VbaunFNZ3wMeE5EgpwZ0Ry2b1fRcB/B06L8mnhMHfEVkJPAuni/Oz6rZ7O/AbJrm7Kg1eL6c13itW+usi6+hhlObd4CrnQ7oAOBxvvo/vwj4ldOnFQ48guc9P+sL4Id82Vy1qspyfdTnM1HVa8AsEblJRPzEc+LE2eatk8CA82xb2z6aRrAEYr7CaSd+GXjE6UCeB/wPnl/7J4D/wvO58cHz5ZaKpz9hOnC/8xzLgTfxnJm0FVhynpf8HM9ZR+kiUt9fteD5QgvG0/zxCp4vjJIGPM/Z53oOzxfMGSART5PHdapaWbWwquao6gptmkl1vsDTeex9BtkaZ12DEpTTJ/UD4HU8tZFTQLJXkSfw9GklALuAbc4675i6eb1+1eX6xFKfz0TVbY/j6Xf7BZ7P2g48HeIAzwMjnCbW6vqhattH0whiE0qZ9kRE/gj0UtU6nY1Vy3P9BrgeuERVcxsdnDHtjNVATJsmIsNEZIxzWu4kPKdivt8Uz62qjwIL8PSBGGOqsBqIadNEZCKeZqveeNrDFwBPNlGzkjHmPCyBGGOMaRBrwjLGGNMgfrUXaT/Cw8O1X79+bodhjDFtytatW7NUtepoDB0rgfTr14/4+Hi3wzDGmDZFRI5Vt96asIwxxjSIJRBjjDENYgnEGGNMg1gCMcYY0yCWQIwxxjSIJRBjjDENYgnEGGNMg1gCaSGf7krjUGZB7QWNMcZRWaks2nycnMJSt0OpliWQesg8XcI/ViRRWFJee2Ev6w5mcf9r23hscXNNGW6MaY9WJ2Xy0Hu7eODtndR33MLP9pxk2e702gs2giWQevj1B4n8ZfkB/rXyYJ23KSgp55fvJOAjsCYpi+PZZ5oxQmNMe7Jo83F8BD7fl8G721LqvF1KbhE/XLSNH72+vVlbPiyB1NHn+06ydHc6Ed068dyaI3VOBH/4ZC+peUX845YJ+Ai8seV4M0dqjGkPMvKL+WxvBndP68/EfqH85qPdpOcV12nbJz/dhyp08vfh1x8k1rv2UleWQOqgqLSCRz7czaDIrrx3/1T8fIXffbLnvNsUl1Xwt8+SeG3Tce6+qD9XjYnismE9eSs+mbKKr82OiqqycO0RLvvLKo5mFTbXrhhjWomH39/Fbc9tJON09Unh7a3JVFQqt0yK5akbx1JWUcm9r8SzP/30eZ9385EcPtqZyvemD+SXc4ex/lA2i3emNscuWAKpi398nkTyqSKeuG4UfcKC+MGlg1i2+ySrD2R+rWxeURkfbE9h7jOrefqzA1w1JooH5gwF4NbJfcgqKGHF3pNf2aasopL/eT+Rx5fs4XBmIb/+8Ku/GMqrSTjGmLbF+/941f4MXtt0nHUHs7nun+vYk5r/lbJnO8+nDOjBgIiu9A/vwtM3jeN4zhmu/PsaHv9oD/vS879WsyirqOTxJbuJCg7k/ukDuXVSLGNjgvntkj3kFZU1+T65OhqviMwF/gb4As+p6pNVHu8EvAxcAGQD31LVo85jD+GZvrQC+LGqLmuuOCsqlZviYrhwQA8A7p7Wn3e2JvP917ax4I4LmDownMSUPJ5atp/1B7Mor1QGhHfh5e9M4pIhX46APH1IJFHBgfzj84McyiykslLZd/I0246dIi2vmO/PGEh41048vmQPSxLSmD40gh8v2s7BjAL+87NLCAqo/e3al57PvrTTXDc+urkOhzEd3omcM6xJyuLWybF1Kv/vVQd5dtUh/nbLeKYM6MEjH+5mQEQX/vzNsXz/1W1c9+91jOsTwtiYYEK7BJCRX0LyqSJ+OXfYuee4YnQUkwf04E/
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADt0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjByYzMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy9h23ruAAAgAElEQVR4nOy9d5wcd33//3pvL7d7e10n3alasmS5YFs2NjYY3DAYMNWhBAgJIbSQQvmFEAIJIZhOiGkGvgEDMRhjSmyDK7Zxt2TZlmVJVtcVXW/b6+f3x2c+s5+Zndlyt3t70n2ej8c97m52duazszOfd39/iDEGhUKhUChkHM0egEKhUCiWHko4KBQKhaIEJRwUCoVCUYISDgqFQqEoQQkHhUKhUJSghINCoVAoSlDCQWELEd1PRO8t8/p3iejTizmm5Q4R/Z6I3l3H45X9jhXLFyUcThCI6AgRJYkoRkQjRPQjImpZxPP/BRE9JG9jjL2fMfa5BpxrLREx7bPGtM/+T9LrREQfJ6L92jU5RkT/SUQeaZ8face4xnTsr2vb/0L6XHnpXDEiut5mXPdr7z3LtP3X2vaX1/M6WMEYexVj7MfS2B+q9J5GUOk7Upz4KOFwYvFaxlgLgBcBOBvAJ5s8nkYT0T7v2wD8KxFdpW3/JoD3AXgXgBCAVwG4HMDPTe9/QdsHAEBELgDXAjho2u9RxliL9PPhMmMyH7MDwIUAxmv9cCcJ4jt6M4BPE9EVzR6Qoj4o4XACwhgbAXAnuJAAABDRBUT0CBHNENEzsharaZiHiChKRIeJ6B3a9s8S0U+l/YQ26JLPR0RbAHwXwIWaljijbf8REf2H9vfLiWiQiD5KRGNEdJyI3iMdo4OI/o+I5ojoSSL6j2q1XsbYowB2AzidiDYC+CCAdzDGHmWM5RhjuwG8CcDVRHSJ9Nb/A3AxEbVp/18F4FkAI9Wc14afAfgzInJq/78NwK8BZKTPej4RPap9F8eJ6HqTVXMlEe0jolki+jYRPSBcO8IaIKKvENG09n29Snrv/UT03jLficFNZLYuiOgKItqrnft6ACR/OCL6SyLao537TiJaU81FYYxtB/+O5HvS8lia5fd17T6ZI6JdRHS69tqPiLsr79bu1wfkMRDRS7T7Z1b7/RLTtfkcET2svfcuIurUXvMR0U+JaFL7Xp4koh7ttVYi+qH2XQ1p96b4fpctSjicgBBRH7i2fED7fxWA2wH8B4B2AB8D8Csi6iKiILim/SrGWAjASwA8Xcv5GGN7ALwfRQ07YrPrCgCtAFYB+CsA35Im5m8BiGv7vFv7qeazEhFdBGArgJ0ALgMwyBh7wjTGAQCPAbhS2pwC8FsAb9X+fxeAG6s5bxmGATwvncfqmHkA/wCgE9yquAxcoEGbrG4Bt/o6AOwD/05kXqxt7wTwJQA/JCLDJF7Dd6KjnftWAP+iHfsggIuk168B8M8A3gigC8CfANxU6bjaey8AcDqK92S5Y10J4GUANoHfL9cCmJQO9w4An9PG+DS4QAYRtYPf598Ev3ZfA3A7cetN8HYA7wHQDcAD/iwA/H5rBdCvvff9AJLaaz8CkANwCrhFfiWAZR+HUcLhxOI3RBQFMABgDMBntO1/DuAOxtgdjLECY+xuANsBvFp7vQCudfsZY8c1TbsRZAH8O2Msyxi7A0AMwKmaFvYmAJ9hjCUYY88D+HEVx5sAMAXgBwD+iTF2L/iEcdxm/+PgE5HMjQDeRUQRAJcA+I3F+y7QtEnxc0GFcYljbgZ3qzwqv8gY28EYe0yzao4A+J52boB/J7sZY7cyxnLgE53ZkjnKGPs+YywPfp16AfRUGFM1iHPfwhjLAviG6dzvB/AFxtgebWz/CeBFFayHCSJKAngUwLdRvL7ljpUFdwduBkDaPvJ3ejtj7EHGWBrAp8Cto34AVwPYzxj7iXZtbwKwF8Brpff+D2PsBcZYEsDNKFoyWXChcApjLK99R3Oa9fBqAH/PGIszxsYAfB1FhWLZooTDicXrNe3/5eAPVqe2fQ2At8gTHICLAfQyxuIA/gz8YT1ORLdrk1ojmNQmAkECQAv4hO0CF2oC+W87OhljbYyxLYyxb2rbJsAnSyt6tdd1GGMPaef/FIDbtEnDzGOMsYj081iFcd0K4FIAHwbwE/OLRLSJiG4jnjgwBz4xiu9qJaTPznjny0HTIUak1xPan/VIPrA6t/w9rAHwX9I9NAXudlpV5pid2tg+Cn5fuisdizF2H4Drwa3JMSK6gYjC0jHlMca0967Ufo6azn/UND5Z2In7D+Df050Afk5Ew0T0JSJya+N0gz8bYqzfA7c8ljVKOJyAMMYeADeFv6JtGgDwE9MEF2SMXaftfydj7ArwyXMvgO9r74sDCEiHXlHutAsY8ji42d4nbeuf57HuA9BPROfLGzXN8gIA91u856fgk9dCXUoA9An79wA+AAvhAOA74Nd5I2MsDO5eEW6h45Cug+Yu6is5QpVDsdhW7js9Dum6a+eWv4cBAH9juo/8jLFHyg6Ca+JfA3fjfbCaYzHGvskYOxfAaeDupY9Lh5TH2ALuKh3WfsxWzGoAQ+XGp50vyxj7N8bYaeBuvNeAuwQHAKTBFRExzjBjbGulY57sKOFw4vINAFcQT6v8KYDXEtEricipBd9eTkR9RNRDRNdosYc0uKunoB3jaQAvI6LVRNSK8tlPowD6SAqsVovmHrkVwGeJKKBZLu+q8Da7Y70AHoj9GfEgvJOItgL4FYBHANxj8bZvArgCwIPzOacN/wzgEs1tZCYEYA5ATPusH5Beux3AGUT0euKB/w+hvFAuh9V38jSAN2rX+RTw2I987q1E9Ebt3B8xnfu7AD6pXU8RqH1LDeO5DsAniMhX7lhEdB4RvVjT3OPgQqUgHefVRHSx9rk+B27ZDQC4A8AmIno7EbmI6M/AhcttlQZGRK8gojM0F+ccuJupoLmz7gLwVSIKE5GDiDaQMbFhWaKEwwkKY2wcXBP+V+3BEQHAcXBt6OPg368DwD+Ca11T4L7vD2jHuBvAL8AzeHag/EN2H3g2yggRTZTZz44PgwcER8C17ZvAhdV8+DB4HOKn4K6D58DdC69njBXMOzPGphhj92pulLrAGBvWXFZWfAw8MBoFt9J+Ib1vAsBbwAPNk+CT23bM71pYfSdfB8+cGgWPV/zM4tzXaefeCOBh6fVfA/giuOtlDvy66plSVXA7gGkAf13hWGHw6zIN/r1NAviydJz/BY+nTQE4FzymBsbYJLjG/1HtPZ8A8Brtc1ViBXgiwByAPQAeQNHqexd48Pp5bUy3wN51uWygOj4vCkXVENEXAaxgjC242peI/g3AGwC8jDE2s+DBLSJE5ACPObyDMfbHZo+n2RDRj8Cz0f6l2WNZ7ijLQbEoENFmIjpTS009H9zd8et6HJsx9hkAN4DHHJY8mvsvQkReFOMRlYLgCsWi4qq8i0JRF0LgrqSV4C6Pr4LXINQFxphly4slyoXgrhPhyni9TRaVQtE0lFtJoVAoFCUot5JCoVAoSmiaW0nLS78RvPKTAbiBMfZfWon8LwCsBXAEwLWMselyx+rs7GRr165t6HgVCoXiZGPHjh0TjDFzVwEATXQrEVEveAXvU0QUAk+lfD2AvwAwxRi7jngL4DbG2P9X7ljbtm1j27dvb/iYFQqF4mSCiHYwxrZZvdY0t5LW4+cp7e8oeO7xKvB8fdF358fgAkOhUCgUi8iSiDkQ0VrwboiPA+iRmnCNwKbhGBG9j4i2E9H28fHl2kpfoVAoGkPThYPWO+VX4F0R5+TXtIpWS78XY+wGxtg2xti2ri5Ll5lCoVAo5klThYPWW+VXAH7GGLtV2zyqxSNEXGKsWeNTKBSK5UrThIPWEfKHAPZoHR0Fv0NxIZh3o46FUgqFQqGojmZWSF8E4J0AdhGRWJnsn8Gbgt1MRH8F3pTr2iaNT6FQKJYtTRMOWkdLsnn5ssUci0KhUCiMND0g3UwOjMXwtbtfwEP7J6DaiCgUCkWRZd14b8/xOVx/334UGHD1mb34+rUvgse
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.optimize as opt\n",
"import scipy.signal as sig\n",
"\n",
"ntaps = 64\n",
"N = 4\n",
"\n",
"#optimize for 16 filter coefficients:\n",
"xmin = opt.minimize(optimfuncQMF, ntaps*np.ones(ntaps), method='SLSQP', tol=1e-8)\n",
"xmin = xmin[\"x\"]\n",
"\n",
"err = optimfuncQMF(xmin)\n",
"print(err)\n",
"\n",
"#Restore symmetric upper half of window:\n",
"h = np.concatenate((xmin, np.flipud(xmin)))\n",
"plt.plot(h)\n",
"plt.title('Resulting PQMF Window Function')\n",
"plt.xlabel('Sample')\n",
"plt.ylabel('Value')\n",
"plt.show()\n",
"\n",
"f, H = sig.freqz(h)\n",
"plt.plot(f, 20*np.log10(np.abs(H)))\n",
"plt.title('Resulting PQMF Magnitude Response')\n",
"plt.xlabel('Normalized Frequency')\n",
"plt.ylabel('dB')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAEWCAYAAADy9UlpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADt0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjByYzMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy9h23ruAAAgAElEQVR4nO3deZwcZZnA8d/Tc2Yyk5kkM5lJJjfJ5CCBQEIAQcADBeTwQA5RjOKyuqJ47Srigcequ+t6gAoiIiAIciorUTmU+0xCIPdMCAk5JpPJNWfmfvaP9+2k05nuubqnumee7+czn+muqq56qrq6nqr3festUVWMMcaY4SAUdADGGGPMYLGkZ4wxZtiwpGeMMWbYsKRnjDFm2LCkZ4wxZtiwpGeMMWbYSHrSE5GbROSbyV5OMoiIisgM/zrueojI10XklsGLLrWIyAgR+T8RqROR+4KOJ5FEZLOIvDvoOGJJ9G9MRG4Tke8nan4R8y0VkadFpEFE/jfR8+9HPGtE5IxeThtzHxCRM0RkWwLjukxEHk3U/BJFRJ4UkU/FGDfVHy8zByGOAS2rx6QXeeCPGHadiNzZmwWo6qdV9Xv+cwPeOUSkQkTuE5Hd/gD7uoh8SUQyBjLfnvS0Hqr6A1XtdodIJBGZKCIPRKz/ahFZkuzl9sKFQCkwVlU/HD1SRIpE5FYR2ekPepUi8rXBDzOxfIJQEbkgavhP/fAlyY4h0b+xJLoS2A2MUtUvR48UkX/3+3ODiLwpIv8ea0YRB76lUcPvFJHrehOMqh6tqk/2bRWST1XvUtX3JGPeInKqiDzvjx17ReQ5ETkhGctKVWlVvCkiRwEvAVuB+apaCHwYWAQUBBnbIPo9bv2nAGOBjwE1gUbkTAEqVbUjxvifAvnAHKAQOB/YOEixHZSkM9FK4PKoZVwEvJGEZaWzKcBajd0jhuC242jgLOAqEbmkh3meKCJvS2CMgUrmlZKIjAL+AtwAjAHKge8ArclaZkpS1bh/gAIzooZdB9zpX58BbAO+DOwCqoFPREx7G/B9YCRwAOgCGv3fBKAZd3UQnv54oBbI6iaWO4FHeoj3fGANsB94EpgTMW4z8BXgdaAO+COQGzH+3338O4BPRq57L9bj4DYZaBw9rF8jsCDGuDOAbVHDNgPvjvje7vPbsQFYBVQA1/jvbivwnjjLnuPXZb9ft/P98O8AbUC7j++Kbj67Gnh/nHmfCaz32+MXwFPAp6L3N/9+qv9uMv37TwDr/DptAv41epsAXwV24k4aQsDXcElpD3AvMCbiMx8Dtvhx10Zuw27ivg34Me7EY7Qfdi7wV+BZYIkfdhTwDz/P3cBdQFHUfv+qX4f7/D7x/QT9xm4Lz6u7/QQ4Dljhl/1H4J6o6c8FVvrv/XngmDjf49uAV/z3+ArwtogY2/1+0hhre0bN63rghhjjwvvAV4F/Rh0jrutN7Bz+2xgB3A7s8/vSf0Rto83E+M1GfD9f99/tZuCyiM8WAnfgjmtbgG8AIT9uCfAc7qRwj/8elwDPRh2DPw1U+fX4JSB+XAbwv365bwJXEfHbiNpmi4D9cbb3dcT/nT0J/BB4GagH/oz/3URMeyXu+FkNfCViXouBF3z81bjfeHYf1vHHfh03AZ+NimuJH97gt8FlsdZRVRN2pVeG+2LLgSuAX4rI6MgJVLUJOBvYoar5/m8HbkNeFDHpx4B7VLW9m+W8G7g/VhAiUgHcDXwBKAGWAv8nItkRk12EO4ucBhyD22CIyFm4nfpMYKZf1hHirEdC4uiFF3Hb9xIRmdzLz0Q6D3fgH407yP4dlwTKge8Cv+7uQyKSBfwf8CgwDvgccJeIzFLVbwM/AP7ot8dvY8T9nyLyCRGZGTXvYuBB3MGgGJeMTunDOu3CHdxG4RLgT0Xk+IjxZbgz2ym4H+XngPcDp+OSwj7cjwwRmQvciNsPJ+Cupif2sPwW3AEgfFVyOe4gd9hq4g4YE3AnD5NwBxn8fvEQLjGMwe07H4j6/EB+YzH5Zf8Jt0+MwSXcD0WMPw64FfhX3Lb4NfCwiOR0M68xwCO4ZDUW+AnwiIiMVdUluET/3z6ux3uIS4C3406u4vkVUNFdfVtfYge+jTtwT8cdAz7azTTxfrNluH23HPg4cLOIzPLjbsB9d9Nx+9zluP007ETcQbsU+M8Y63kucIJf7kXAe/3wf8F95wtwJ07vj/F5cCUSnSJyu4icHb3/9NLluAuC8UAH7ruO9A7c8fM9wFcjvpdO4Iu4bXQy8C7g36I+G28dz8WdnC3CVaUAICIjfQxnq2oB7qRrZbwVSFTSawe+q6rtqroUdyY3q4fPhN2O38F8vdyluB9gd8bizhJiuRh3JfiYT5o/xp3B/UlEduEOONer6g5V3Ys7iC/wn70I+J2qrvYHj+viLOcJYJyIrBSRh/sQR2QxTKw4evJh4Bngm8CbPoa+lMk/o6p/V1cMeR8uKf/Ix3kPMFVEirr53Em44skfqWqbqv4DV1RyaS+X+zncQe8qYK2IbBSRs/24c4A1qnq/j+NnuKuyXlHVR1T1DXWewiXmt0dM0gV8W1VbVfUA7ozyWlXdpqqtuO/6Ql+0dCHwF1V92o/7pv98T+4ALvfb7nRcIomMcaPfH1pVtRaXEE73o08CMnH7RLuqPog7m440kN9YPCcBWcDP/Lzvx12hhV0J/FpVX1LVTlW9HVccdlI383ofUKWqv1fVDlW9G3f1fl4/4roOd3z6XQ/THcAliu4a3vQl9ouAH6jqPlXdxpEHc+j5N/tN//0+hUv+F/lj2iXANaraoKqbcVdmH4v43A5VvcFvswMx1vNHqrpfVd8C/snhx62f+315H/CjGJ9HVeuBU3FXSb8BakXkYREpjfWZbvw+4hj5zYh1DPuOqjap6ircd3epX/ZyVX3Rr+Nm3AnI6VHzjreOP1PVrX7b/zDqc13APBEZoarVqhr3RKk3Sa8T96OIlIX7EYbt0cPrcppxB8je+DMwV0Sm4c6w6lQ1+gd/cDm4M4xYJuCKDwBQ1S5ckd0y3BkaHH4wjYxzgp82bAuxtQG7VHWBqp7fhzjKI6aJFUdc/kf5NVU9GndmuBKX1KU3n+fw+r8DwG5V7Yx4T4xYJgBb/bqEbeHwdYoX9wF1jX0W4k5e7gXu81cHh217dWUWW7uf05H8WeuLvmJ+Py6JFkdMUquqLRHvpwAPich+P/063H5e2k0sTbj9rqf1exZ3AnEtLmkedvAS13LxHhHZLiL1uGK4cIwTgO1+vcOi138gv7F4ult25L4/BfhyeFv57TXJf667eUX/bnq9j4SJyFW4K4r3+ROPntwClIpIdHLta+yR27y7/S/eb3af31fCtvh5FuOOl1uixkVuk97s6709bsWdl6quU9UlqjoRmOc//7NeLL+7+W/BrVtxnPET4GADxL/4hmz1uJKhyM9BP47NfptfjDuRrRaRR0RkdrwV6E3Sewt32R9pGvGTQixHVGD7g9G9uKu9jxH7Kg/gcSKKXrqxA7ejAweLSCbh6lL2Rk7oG8V8FDhbRJ7BbeRJEZPEKzrs6dEUseLY3sPn+kRVd+OuIifgiqaagLyI5WbgDsSJsAOYJCKR+8xk+rFO/ozzB7g6qGm4q/eD2z5ie4Udtl64oqTwtDnAA7jtUKqqRbji5MiTgOjvayuuOKQo4i9XVbd3E0seLkn3xp24erfook1w66u4BlijcPteOMZqoDzqxGUS/dPdvhlz+8VYduS+vxX4z6htleev4qIdtt9HzKvX+4iIfBJX3/ouf8XVI1Vtw9Urf4/Dv/e+xF7N4cXYfd3+o31RW9hk3PbYjbtAmBI1LnKbDORRN/2OW1XX44rU5/lB8faT7uY/Gbduu+OMDxev34i76p/p9/+vc/h3Fc9hv0mijs2+5OpM3AXRetxVbEy9SXp/BL4hrql8yJfRnkecurU4aoCxIlIYNfwOXPn4+cRPet8G3iYi/yMiZQAiMsM3Uy7CJc/3ici7fB3Ul3HFGc93M6+bcQfHv+Lq8qYDS0Rkrj/QfTtOHNnABBF5RUS6K0PvSxx
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"N = 4\n",
"f, H_im = sig.freqz(h)\n",
"posfreq = np.square(H[0:512//N])\n",
"negfreq = np.flipud(np.square(H[0:512//N]))\n",
"plt.plot((np.abs(posfreq) + np.abs(negfreq)))\n",
"plt.xlabel('Frequency (512 is Nyquist)')\n",
"plt.ylabel('Magnitude')\n",
"plt.title('Unity Condition, Sum of Squared Magnitude of 2 Neighboring Subbands')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"b = sig.firwin(80, 0.5, window=('kaiser', 8))"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAEWCAYAAADy9UlpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADt0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjByYzMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy9h23ruAAAgAElEQVR4nO3deZwcZZnA8d/Tc2Yyk5kkM5lJJjfJ5CCBQEIAQcADBeTwQA5RjOKyuqJ47Srigcequ+t6gAoiIiAIciorUTmU+0xCIPdMCAk5JpPJNWfmfvaP9+2k05nuubqnumee7+czn+muqq56qrq6nqr3festUVWMMcaY4SAUdADGGGPMYLGkZ4wxZtiwpGeMMWbYsKRnjDFm2LCkZ4wxZtiwpGeMMWbYSHrSE5GbROSbyV5OMoiIisgM/zrueojI10XklsGLLrWIyAgR+T8RqROR+4KOJ5FEZLOIvDvoOGJJ9G9MRG4Tke8nan4R8y0VkadFpEFE/jfR8+9HPGtE5IxeThtzHxCRM0RkWwLjukxEHk3U/BJFRJ4UkU/FGDfVHy8zByGOAS2rx6QXeeCPGHadiNzZmwWo6qdV9Xv+cwPeOUSkQkTuE5Hd/gD7uoh8SUQyBjLfnvS0Hqr6A1XtdodIJBGZKCIPRKz/ahFZkuzl9sKFQCkwVlU/HD1SRIpE5FYR2ekPepUi8rXBDzOxfIJQEbkgavhP/fAlyY4h0b+xJLoS2A2MUtUvR48UkX/3+3ODiLwpIv8ea0YRB76lUcPvFJHrehOMqh6tqk/2bRWST1XvUtX3JGPeInKqiDzvjx17ReQ5ETkhGctKVWlVvCkiRwEvAVuB+apaCHwYWAQUBBnbIPo9bv2nAGOBjwE1gUbkTAEqVbUjxvifAvnAHKAQOB/YOEixHZSkM9FK4PKoZVwEvJGEZaWzKcBajd0jhuC242jgLOAqEbmkh3meKCJvS2CMgUrmlZKIjAL+AtwAjAHKge8ArclaZkpS1bh/gAIzooZdB9zpX58BbAO+DOwCqoFPREx7G/B9YCRwAOgCGv3fBKAZd3UQnv54oBbI6iaWO4FHeoj3fGANsB94EpgTMW4z8BXgdaAO+COQGzH+3338O4BPRq57L9bj4DYZaBw9rF8jsCDGuDOAbVHDNgPvjvje7vPbsQFYBVQA1/jvbivwnjjLnuPXZb9ft/P98O8AbUC7j++Kbj67Gnh/nHmfCaz32+MXwFPAp6L3N/9+qv9uMv37TwDr/DptAv41epsAXwV24k4aQsDXcElpD3AvMCbiMx8Dtvhx10Zuw27ivg34Me7EY7Qfdi7wV+BZYIkfdhTwDz/P3cBdQFHUfv+qX4f7/D7x/QT9xm4Lz6u7/QQ4Dljhl/1H4J6o6c8FVvrv/XngmDjf49uAV/z3+ArwtogY2/1+0hhre0bN63rghhjjwvvAV4F/Rh0jrutN7Bz+2xgB3A7s8/vSf0Rto83E+M1GfD9f99/tZuCyiM8WAnfgjmtbgG8AIT9uCfAc7qRwj/8elwDPRh2DPw1U+fX4JSB+XAbwv365bwJXEfHbiNpmi4D9cbb3dcT/nT0J/BB4GagH/oz/3URMeyXu+FkNfCViXouBF3z81bjfeHYf1vHHfh03AZ+NimuJH97gt8FlsdZRVRN2pVeG+2LLgSuAX4rI6MgJVLUJOBvYoar5/m8HbkNeFDHpx4B7VLW9m+W8G7g/VhAiUgHcDXwBKAGWAv8nItkRk12EO4ucBhyD22CIyFm4nfpMYKZf1hHirEdC4uiFF3Hb9xIRmdzLz0Q6D3fgH407yP4dlwTKge8Cv+7uQyKSBfwf8CgwDvgccJeIzFLVbwM/AP7ot8dvY8T9nyLyCRGZGTXvYuBB3MGgGJeMTunDOu3CHdxG4RLgT0Xk+IjxZbgz2ym4H+XngPcDp+OSwj7cjwwRmQvciNsPJ+Cupif2sPwW3AEgfFVyOe4gd9hq4g4YE3AnD5NwBxn8fvEQLjGMwe07H4j6/EB+YzH5Zf8Jt0+MwSXcD0WMPw64FfhX3Lb4NfCwiOR0M68xwCO4ZDUW+AnwiIiMVdUluET/3z6ux3uIS4C3406u4vkVUNFdfVtfYge+jTtwT8cdAz7azTTxfrNluH23HPg4cLOIzPLjbsB9d9Nx+9zluP007ETcQbsU+M8Y63kucIJf7kXAe/3wf8F95wtwJ07vj/F5cCUSnSJyu4icHb3/9NLluAuC8UAH7ruO9A7c8fM9wFcjvpdO4Iu4bXQy8C7g36I+G28dz8WdnC3CVaUAICIjfQxnq2oB7qRrZbwVSFTSawe+q6rtqroUdyY3q4fPhN2O38F8vdyluB9gd8bizhJiuRh3JfiYT5o/xp3B/UlEduEOONer6g5V3Ys7iC/wn70I+J2qrvYHj+viLOcJYJyIrBSRh/sQR2QxTKw4evJh4Bngm8CbPoa+lMk/o6p/V1cMeR8uKf/Ix3kPMFVEirr53Em44skfqWqbqv4DV1RyaS+X+zncQe8qYK2IbBSRs/24c4A1qnq/j+NnuKuyXlHVR1T1DXWewiXmt0dM0gV8W1VbVfUA7ozyWlXdpqqtuO/6Ql+0dCHwF1V92o/7pv98T+4ALvfb7nRcIomMcaPfH1pVtRaXEE73o08CMnH7RLuqPog7m440kN9YPCcBWcDP/Lzvx12hhV0J/FpVX1LVTlW9HVccdlI383ofUKWqv1fVDlW9G3f1fl4/4roOd3z6XQ/THcAliu4a3vQl9ouAH6jqPlXdxpEHc+j5N/tN//0+hUv+F/lj2iXANaraoKqbcVdmH4v43A5VvcFvswMx1vNHqrpfVd8C/snhx62f+315H/CjGJ9HVeuBU3FXSb8BakXkYREpjfWZbvw+4hj5zYh1DPuOqjap6ircd3epX/ZyVX3Rr+Nm3AnI6VHzjreOP1PVrX7b/zDqc13APBEZoarVqhr3RKk3Sa8T96OIlIX7EYbt0cPrcppxB8je+DMwV0Sm4c6w6lQ1+gd/cDm4M4xYJuCKDwBQ1S5ckd0y3BkaHH4wjYxzgp82bAuxtQG7VHWBqp7fhzjKI6aJFUdc/kf5NVU9GndmuBKX1KU3n+fw+r8DwG5V7Yx4T4xYJgBb/bqEbeHwdYoX9wF1jX0W4k5e7gXu81cHh217dWUWW7uf05H8WeuLvmJ+Py6JFkdMUquqLRHvpwAPich+P/063H5e2k0sTbj9rqf1exZ3AnEtLmkedvAS13LxHhHZLiL1uGK4cIwTgO1+vcOi138gv7F4ult25L4/BfhyeFv57TXJf667eUX/bnq9j4SJyFW4K4r3+ROPntwClIpIdHLta+yR27y7/S/eb3af31fCtvh5FuOOl1uixkVuk97s6709bsWdl6quU9UlqjoRmOc//7NeLL+7+W/BrVtxnPET4GADxL/4hmz1uJKhyM9BP47NfptfjDuRrRaRR0RkdrwV6E3Sewt32R9pGvGTQixHVGD7g9G9uKu9jxH7Kg/gcSKKXrqxA7ejAweLSCbh6lL2Rk7oG8V8FDhbRJ7BbeRJEZPEKzrs6dEUseLY3sPn+kRVd+OuIifgiqaagLyI5WbgDsSJsAOYJCKR+8xk+rFO/ozzB7g6qGm4q/eD2z5ie4Udtl64oqTwtDnAA7jtUKqqRbji5MiTgOjvayuuOKQo4i9XVbd3E0seLkn3xp24erfook1w66u4BlijcPteOMZqoDzqxGUS/dPdvhlz+8VYduS+vxX4z6htleev4qIdtt9HzKvX+4iIfBJX3/ouf8XVI1Vtw9Urf4/Dv/e+xF7N4cXYfd3+o31RW9hk3PbYjbtAmBI1LnKbDORRN/2OW1XX44rU5/lB8faT7uY/Gbduu+OMDxev34i76p/p9/+vc/h3Fc9hv0mijs2+5OpM3AXRetxVbEy9SXp/BL4hrql8yJfRnkecurU4aoCxIlIYNfwOXPn4+cRPet8G3iYi/yMiZQAiMsM3Uy7CJc/3ici7fB3Ul3HFGc93M6+bcQfHv+Lq8qYDS0Rkrj/QfTtOHNnABBF5RUS6K0PvSxx
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"f, H_im = sig.freqz(h)\n",
"posfreq = np.square(H[0:512//N])\n",
"negfreq = np.flipud(np.square(H[0:512//N]))\n",
"plt.plot((np.abs(posfreq) + np.abs(negfreq)))\n",
"plt.xlabel('Frequency (512 is Nyquist)')\n",
"plt.ylabel('Magnitude')\n",
"plt.title('Unity Condition, Sum of Squared Magnitude of 2 Neighboring Subbands')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"text/plain": [
"(63,)"
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b.shape"
]
},
{
"cell_type": "code",
"execution_count": 102,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEWCAYAAACqitpwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADt0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjByYzMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy9h23ruAAAgAElEQVR4nO3deXycZbn4/8+VTPakSZomXdOk0JRSCrSQtqAIiKBFloIilIMCfjlf9Kscj4eDCi6oHDgH1N9BPeKCiIKsPSxSZKkgsimUpjultE3bpEmXrE2bpVkmc/3+eJ4p0yF7ZuaZSa/36zWvzNzPdt/JZK65l+e+RVUxxhhjIinJ6wwYY4wZeyy4GGOMiTgLLsYYYyLOgosxxpiIs+BijDEm4iy4GGOMiTgLLiZuiMirIvLPA2z/tYh8L5Z5ShQiMl1E2kQkeYTHV4nIuZHOVzwRkU0icrbX+ThaWHAxfXI/bA65H1j7ROQPIpIdw+tfKyJvhqap6pdV9T+icK1SEVG3rG1u2W8O2S4i8g0R2eb+TnaJyH+KSGrIPn9wz7Ek7Nx3u+nXhpSrN+RabSLyiz7ydKWIbA5Le6mftJtVdZeqZqtqb4R+LRElIj8QkZ6wcn8zitf7g4jcHpqmqieo6qvRuqY5kgUXM5CLVDUbmAfMB27xOD/RlueW90rgVhFZ7Kb/HLgeuBrIAc4HzgUeCzt+q7sPACLiAy4Htoft95YbCIKPG/rIy+vAbBEpDDnXyUBGWNrp7r6J4PGwcv/I6wyZ6LHgYgalqvuAFThBBgAROU1E/iEiLSKyPrS5wf12vkNEWkVkp4hc5ab/QEQeCtkvWGPwhV5PRI4Hfg2c7n7DbXHTD38bFZGzRaRWRP5dROpFZK+IfDHkHAUi8qyIHBSRVSJye3hNaIDyvgVsAuaKSBnwFeAqVX1LVf2qugn4LHCBiJwVcuizwBkiku++XgxsAPYN5bphedgN7ADOdJNOcfP0WlhaErAq/HfpNjH+h4j83f07/EVEJoT8fr4gItUi0iQi3wm9toikichPRWSP+/ipiKS5214Tkc+6zz/qXvMC9/UnRGTdcMo52HtiCOU4I+R9WOO+964HrgK+6b5/nnX3Pdz0N0gZB3xvmaGx4GIGJSLTcL6tV7qvpwLPAbcD44GbgCdFpFBEsnC+6Z+vqjnAR4BhfeCo6mbgy3zwDT+vn10nAbnAVOA64J6QD/Z7gHZ3n2vcx1DKKiLyUeAEYC3wCaBWVd8Jy2MN8DbwyZDkTuAZYKn7+mrgwaFctx+v80EgORN4A3gzLO1tVe3p5/h/Ar4IFAGpOH8nRGQO8CvgC8AUoACYFnLcd4DTcL5MnAwsBL7rbnsNONt9fhZHBsCz3O2R1l85SoAXgP8BCt38rlPVe4GHgR+575+L+jjnQGWEgd9bZggsuJiB/ElEWoEaoB74vpv+eeB5VX1eVQOq+hJQAXza3R7A+dafoap73W/60dAD3KaqPar6PNAGHCdOp/Znge+raoeqvgc8MITzNQLNwH3Azar6V2ACsLef/ffifKiFehC4WkTycD5s/9THcae537SDj9P6OX9oLeVjOMHljbC0gT7Mf6+qW1X1ELCMD2qelwF/VtXXVbUL+B7O3yzoKpzfa72qNgA/xAlEwTwFa2tnAv8V8nqw4HJ5WLmnDLDvUMrxT8DLqvqo+x5oUtWhfpEZqIzQz3triOc2WHAxA7vErX2cDczG+aAFKAE+F/pBAZwBTFbVduAKnJrHXhF5TkRmRyl/TarqD3ndAWTjfOD7cIJiUOjz/kxQ1XxVPV5Vf+6mNQKT+9l/srv9MFV9073+d3A+wA/1cdzbqpoX8ni7n/O/DpzkfmM+Dacm9z4w2U07g4H7W0Kb44K/G3BqK4d/H+7frClk3ylAdcjrajcN4C1glohMxPmQfxAodpuqFg6Sn2Vh5d4zwL5DKUcxH+7PGqqBygj9v7fMEFlwMYNS1deAPwA/cZNqgD+GfVBkqeqd7v4rVPU8nA/f94Hfuse1A5khp5400GVHkeUGwM+RTT3FIzzXKzgfngtDE0WkGOcD/9U+jnkI+HdG1ySGqu4A9uAMJtilqm3uprfctGycprnh2kvI70NEMnGaxoL24HyBCJrupqGqHcBq4F+Bd1W1G/gHcCOwXVWPCLZDMJz3RLga4Nh+tg32/um3jCYyLLiYofopcJ6InIzz4XmRiHxKRJJFJN3tBJ0mIhNFZInb99KF05wQbHJZB5wpzj0ZuQw8+qwOmCYhw32Hyh2O+xTwAxHJdGtOVw9yWH/n2oozuOBhcQYxJIvICcCTOB+qL/dx2M+B84jMKK43cD643whJe9NNq+inZjSYJ4AL3c7wVOA2jvwseBT4rtuHNgG4FedvHvQacAMfNIG9GvZ6OIbzngj3MHCuiFwuIj5xBnEEm8zqgGMGOHawMppRsuBihsRtl34QuNXtzF4CfBunllADfAPn/ZSE88G3B6f/4izg/7nneAl4HGcE1WrgzwNc8hWc0VH7RGS434bB+bDLxWlS+SPOh0nXCM4TPNd9OB8+HcC7OM0ol6hqIHxnVW1W1b9qZBZLeg2nIzt0pNsbbtqIgpfbB/ZV4BGcWsx+oDZkl9tx+tA2ABuBNW5aaJ5yQq4f/no4eRnOeyL82F04/Xz/jvNeW4fTOQ/wO2CO22zbV7/XYGU0oyS2WJg5GojIXcAkVR3SqLFBzvVD4FLgTFVtGXXmjBmDrOZixiQRmS0iJ7lDixfiDCd9OhLnVtXvA/fi9LkYY/pgNRczJonIApymsCk47e/3AndGqKnKGDMICy7GGGMizprFjDHGRJxv8F3GvgkTJmhpaanX2TDGmISyevXqRlUNn6UCsOACQGlpKRUVFV5nwxhjEoqIVPe3zZrFjDHGRJwFF2OMMRFnwcUYY0zEWXAxxhgTcRZcjDHGRJynwUVEFovIFhGpFJGb+9h+poisERG/iFwWtq1XRNa5j+Uh6TNEZKV7zsdHMquuMcaY0fEsuLirBd6Ds3zuHOBKd/nVULuAa3Fmbw13SFXnuY+LQ9LvAu5W1Zk4s71eF/HMG2OMGZCXNZeFQKWq7nAXHHoMZxr3w1S1SlU3cOQSrP0SEQHOwVmvApylbS+JXJaNiZ2Obj/LKmro7On1OivGDJuXwWUqRy49W+umDVW6iFSIyNsiEgwgBUBLyPKk/Z5TRK53j69oaGgYbt6Nibpbn9nEN5/YwJ0vvO91VowZtkTu0C9R1XLgn4Cfikh/y532SVXvVdVyVS0vLOxz9gJjPPPiu3t5YnUt08dn8od/VPHGNvsCZBKLl8FlN0euaz7NTRsSVd3t/tyBs8zqfKAJyBOR4LQ2wzqnMfGgvrWTW57ayIlTc3nua2cwsyibm/53PS0d3V5nzZgh8zK4rALK3NFdqcBSYPkgxwAgIvkikuY+nwB8FHjPXavjb0BwZNk1wDMRz7kxUaKqfPOJDXR093L3FSeTk57CT6+YR1NbN9/907vYEhkmUXgWXNx+kRuAFcBmYJmqbhKR20TkYnAWfBKRWuBzwG9EZJN7+PFAhYisxwkmd6rqe+62bwE3ikglTh/M72JXKmNG5+GVu3h1SwPf/vTxzCzKAWDu1Fz+7bxZ/HnDXpav3+NxDo0ZGlssDCgvL1ebFdl4bUdDGxf8/E3KS/N54IsLSUqSw9v8vQGuuPdttta18uLXz2RqXoaHOTXGISKr3b7vD0nkDn1jxgx/b4B/W7aeVF8SP77s5CMCC4AvOYn/vvxkAgHlpmXrCQTsS6GJbxZcjIkDK3c2s76mhe9dOIdJuel97lNSkMU3F8/mrR1NbNx9IMY5NGZ4LLgYEwdW7mgiSeBTJ0wccL/z505y9t/ZFItsGTNiFlyMiQNv72zmhCm55KSnDLhf0bh0ZkzIYuWO5hjlzJiRseBijMc6e3pZV9PCohnjh7T/ohnjeaeqmV7rdzFxzIKLMR5bV9NCtz/AomMKhrT/omPG09rpZ/Peg1HOmTEjZ8HFGI+t3NGMCCwsHWrNxQlCK3da05iJXxZcjPHYyp1NzJ4
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"cutoff = 0.15\n",
"beta = 9\n",
"ntaps = 63\n",
"N = 4\n",
"\n",
"b = sig.firwin(ntaps, cutoff, window=('kaiser', beta))\n",
"w, h = sig.freqz(b)\n",
"\n",
"plt.plot(b)\n",
"plt.title('Resulting PQMF Window Function')\n",
"plt.xlabel('Sample')\n",
"plt.ylabel('Value')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAEWCAYAAADl19mgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADt0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjByYzMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy9h23ruAAAgAElEQVR4nOydd3hVVdb/PyudhN6r9F4SEBTbCIoDCoqAIkrJ1XF8HXVmnFffmZ8z48So4ziOfazojBQLomLFLiAWUBAIvXcIvbeEJOv3xz7n5ubmtoQkN4n78zznufecs84+617C+d6999priapisVgsFkt1IibaDlgsFovFUtZYcbNYLBZLtcOKm8VisViqHVbcLBaLxVLtsOJmsVgslmqHFTeLxWKxVDusuFlCIiIviMi9ZW0b5HoVkQ5BzjURkbkiclREHhORP4vIy865Ns61caW9dwS+Fbl/ed3HYrGUDeX2MLBUfkRkM9AEyAPygZXAFGCiqhYAqOqtkbbnaysiA4BXVbVlGbl7C7APqK1hFmeKyBzn3i+X0b1LdH+LxRJ9bM/NcqWq1gJaAw8DfwL+E12XAtIaWFkRwiIisSW9f3n2GisjQb4ji6XyoKp2+5luwGZgkN+xc4ACoIezPwl40Of8H4FsYCdwM6BAB19bIAU46bRzzNmaO23PAw45bTwDJPi07W3Lz6dJwGkg12lrEHAfpncG0Ma5Ng74O6YXesqxfcax6QJ8ARwA1gCj/dp/HvgYOB7gOwl2/7eBV4EjzndRB/PDIBvY4XwXsU4bscCjmN7fRuB21+dA/xa+n8/Z7w9873x3WcAAn3NzgAeA74CjwOdAQ5/zF/pcuw3wAP2A3a5/jt1IICvI30qx78j5N30H2AtsAn7n93e00PludgOP+/1b3YL5G8oG7va5LhF40jm303mf6JwbAGwH7gL2ONfe6HPtFZjRh6PO9+/b7jBgifMdfA/0ivb/P7uV7xZ1B+wWxX/8AOLmHN8K/MZ5PwlH3IAhwC6gO5DsPNiLiZvzfgCw3a/ds52HdJzzkFsF3OlzPqC4+bft7Hsf/j4PTFco5gA3+9imOA/1G51798aITDeftg8DF2BGM5IivP9p4GrnmhrAu8CLzv0aAz8C/+PY3wqsBloB9YHZRChuQAtgv/PwjgEuc/Yb+XzeDUAnx485wMPOudaYh/31QDzQAEhzzq0ELve557vAXSG+f9/vKBn4CfgbkAC0w4j2YMd+HjDeeV8T6O/3b/WG8z31xIjjIOf8/cB85/trhBGiB3z+pvIcm3jn+zgB1HPOZwMXOe/rAX2c970xYngu5kdGuvN9J0b7/6Ddym+zw5KWQOzEPID9GQ28oqorVPUE5gEcMar6k6rOV9U8Vd2MEYKLz9TZCBgGbFbVV5x7L8b0OK71sXlfVb9T1QJVPRVhu/NU9T0185O1MQ/bO1X1uKruAZ4Axji2o4EnVXWbqh4A/lEC/8cBH6vqx45/X2B6RVf42LyiqmtV9SQwHUhzjt8AfKmqb6jqaVXdr6pLnHOTnbYRkfrAYOD1EH54vyOMKDVS1ftVNVdVNwIv+Xze00AHEWmoqsdUdb5fW5nO97QMeAUjvgBjgftVdY+q7gUygfE+1512zp9W1Y8xPenOPue6iUhtVT2oqouc47cAL6rqD6qar6qTgRzMDy1LNcWKmyUQLTDDd/40x/SAXLYFsAmKiHQSkY9EZJeIHAEeAhqW3s2IaQ2cKyKH3A3zEG3qY1OizxLgmtaY3kS2zz1exPRAoPh3t6UE92kNXOvn/4VAMx+bXT7vT2B6S2B6ihuCtPsqcKWIpGDE9xtVzQ7hh//nbe7n058xAUoAv8L0JFeLyAIRGRairS2Y7wfndUuQcwD7VTXPZ9/3s47CCP4WEflaRM7z8fUuP19b+bVrqWb8rCbBLeERkX4Ycfs2wOlswDf6sVWIpgIFXjwPLAauV9WjInIncE1pfS3BvbcBX6vqZSW4pqT32YbpDTT0e/i6ZFP0+zrL7/xxzFCfi7/wTlXVX5fCx22Y+a9iqOoOEZmHmWsbj/n3CYX/592kqh2DtL0OuF5EYpz23xaRBj4mrTDDtGC+i53O+50YMVoR4Fxo51QXAMNFJB64A9ODbeX4+ndV/Xsk7ViqB7bnZgFARGo7v66nYeZ6lgUwmw7cKCJdRSQZCLWmbTfQQETq+ByrhQkwOCYiXYDflJH7ge7dzmf/I6CTiIwXkXhn6yciXcvqhk6P53PgMee7jBGR9iLiDrtOB34nIi1FpB7w//yaWAKMcXzrS1HRd3tYg0UkVkSSRGSAiESyzOI1YJCIjBaROBFpICJpPuenYIKEegIzSvCRfwSOisifRKSG41cP58cRIjJORBo5Q5iHnGsKfK6/V0SSRaQ7Zi70Tef4G8BfRaSRiDTEzOm9Gs4ZEUkQkbEiUkdVT2P+ztz7vQTcKiLniiFFRIaKSK0SfF5LFcOKm+VDETmK+XX7F+BxzMOmGKr6CfA0JhhiPWbiH0yPxd92NeZBtdEZCmoO3I2ZAzqKeeC86X9dGfEUcI2IHBSRp1X1KPBLzHzQTswQ3j8xkXllyQRMcMVK4CAmmtIdOnwJ+AwT6biI4kJyL9DeuS4Tn7kvVd0GDMcM++3F/Fv9HxH8/1XVrZihurswQ81LgFQfk3cxPaV3nXnUiFDVfMxcZhomUnIf8DImYhRM8NEKETmG+fcY48wHunyN+Rv6CnhUVT93jj+ImU9cCizDfFcPRujWeGCzM+R9K2boGVVdCPwaE5170LmvJ9LPaqmaiKpdj2opHU7PZzkm6izQUJwlCCLSBiMK8dH+7kRkAyaq88sKuFcbKsnntlRvbM/NUiJEZISIJDpDa/8EPrQPqaqLiIzCzKXNirYvFktZYsXNUlL+B7NmaANmsXR5zZtZyhknTdnzwO3O3JjFUm2ww5IWi8ViqXbYnpvFYrFYqh3Vfp1bTEyM1qhRI9puWCwWS5XhxIkTqqpVuvNT7cWtRo0aHD9+PNpuWCwWS5VBRE6Gt6rcVGlltlgsFoslEFbcLBaLxVLtsOJmsVgslmqHFTeLxWKxVDusuFksFoul2lHlxE1EhojIGhFZLyL+mdUtFovFUgFU9mdxlRI3EYkFngUuB7ph6kV1i65XFovF8vOiKjyLq9o6t3OA9U5Je0RkGqYUyMqyvMnJ0ycZ8Y/naJTXm3ZyCSKU21azJtSrV7i1aAF165blp7FYLJYyp0KexWdCVRO3FhQtT78dONffSERuAW4BaNGiBXPmzCnRTfI1n+95iLMSupBWo/SdW1WzheLkSdi7t+ixuDhITIQaNYz41axp9i0Wi6WCiBORhT77E1V1os9+RM/iaFLVxC0inH+EiQApKSk6YMCAErfxP6dv4skfnmT8rd1olNzYK1RlvR07BgcPwoED5nXbNli1CtauhcWL4ZBTw7hLFxg5EiZMgM6dy/DLslgsluLkqWrfaDtxJlQ1cdsBtPLZb+kcK3PS09J5dN6jvL7sde7sfyci5XEXaNIk+LmCAlizBmbNgnffhX/+Ex56CAYPhrvvhkGDyscni8ViCUOFPYtLS5UqeSMiccBa4FLMF7kAuEFVVwS7JiUlRUubW7LfS/04nX+aJbcuKdX1Zc3u3TBxIjz/PGRnG5F79FHo0SPanlksluqEiJxQ1ZQQ50v8LK5oqlS0pFPx+Q7gM2AVML08v8z01HSydmexZFflELcmTeDee2HTJnj8cfjhB0hNhb/+FXJzo+2dxWL5uVDRz+LSUKV6bqXhTHpu+0/sp9ljzbi93+08MeSJMvbszNm/3wxPTpoEffrAjBnQunW0vbJYLFWdcD23qkCV6rlVNA2SG3BV56t4bdlrnM4/HW13itGgAbzyipmP27AB+vaFuXOj7ZXFYrFEHytuYUhPTWfvib18sv6TaLsSlKuvhh9/hIYN4bLL4IMPou2RxWKxRBcrbmEY0mEIjVMaM2nJpGi7EpJOneC77yAtzSwZePPNaHtksVgs0cOKWxjiY+MZ13McH639iH0n9kXbnZDUrw9ffgnnnw/
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"fig, ax1 = plt.subplots()\n",
"ax1.set_title('Digital filter frequency response')\n",
"\n",
"ax1.plot(w / (2 * np.pi), 20 * np.log10(abs(h)), 'b')\n",
"ax1.set_ylabel('Amplitude [dB]', color='b')\n",
"ax1.set_xlabel('Frequency [rad/sample]')\n",
"\n",
"ax2 = ax1.twinx()\n",
"angles = np.unwrap(np.angle(h))\n",
"ax2.plot(w / (2 * np.pi), angles, 'g')\n",
"ax2.set_ylabel('Angle (radians)', color='g')\n",
"ax2.grid()\n",
"ax2.axis('tight')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 105,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"text/plain": [
"(63,)"
]
},
"execution_count": 105,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b.shape"
]
},
{
"cell_type": "code",
"execution_count": 129,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"def optimfuncQMF(x):\n",
" \"\"\"Optimization function for a PQMF Filterbank\n",
" x: coefficients to optimize (first half of prototype h because of symmetry)\n",
" err: resulting total error\n",
" \"\"\"\n",
" N = 2 #4 subbands\n",
" cutoff = 1.5 #1.5\n",
" h = np.append(x, np.flipud(x))\n",
" f, H_im = sig.freqz(h)\n",
" H = np.abs(H_im) #only keeping the real part\n",
" \n",
" posfreq = np.square(H[0:512//N])\n",
" \n",
" #Negative frequencies are symmetric around 0:\n",
" negfreq = np.flipud(np.square(H[0:512//N]))\n",
" \n",
" #Sum of magnitude squared frequency responses should be closed to unity (or N)\n",
" unitycond = np.sum(np.abs(posfreq + negfreq - 2*(N*N)*np.ones(512//N)))//512\n",
" \n",
" #plt.plot(posfreq+negfreq)\n",
" \n",
" #High attenuation after the next subband:\n",
" att = np.sum(np.abs(H[int(cutoff*512//N):]))//512\n",
" \n",
" #Total (weighted) error:\n",
" err = unitycond + 100*att\n",
" return err"
]
},
{
"cell_type": "code",
"execution_count": 131,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3.0\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEWCAYAAABIVsEJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADt0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjByYzMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy9h23ruAAAgAElEQVR4nO3dd3xUVfr48c+TRgglhSQQEkLvHQMIoqCAYEVX17Wsoqvr6vbid7/6dVdd19113aJbv/74KnaxF0SFRQTpJbQQauikkUYSEtLz/P6YC44xIT035Xm/XvPK3Dvnzjz3zmSeOefce46oKsYYY0x9+bgdgDHGmLbJEogxxpgGsQRijDGmQSyBGGOMaRBLIMYYYxrEEogxxpgGsQRiWpyIrBKRe87z+LMi8uuWjKmtEJFYESkQEd8Gbn9URGY1dVytiYjsFpEZbsfREVgC6eCcL5Qi50spXUReFJGuLfj6d4rIWu91qnqfqv62GV6rn4ios68Fzr4/6PW4iMh/iUiSc0yOi8jvRSTAq8yLznPMq/LcTzvr7/Tarwqv1yoQkX9WE9MtIrK3yrrlNax7UFWPq2pXVa1oosPSpETkMREpq7Lfv2zG13tRRJ7wXqeqI1V1VXO9pvmSJRADcI2qdgXGAeOBh1yOp7mFOPt7C/CIiMx11v8duBe4A+gGXAHMAt6osv0BpwwAIuIH3AQcqlJug/Nlf/b2w2piWQ0ME5EIr+caC3Susm6KU7YteLPKfj/ldkCmeVgCMeeoajqwDE8iAUBELhSR9SKSKyI7vZsGnF/Zh0XktIgcEZHbnPWPicirXuXO/vL38349ERkOPAtMcX6p5jrrz/2qFJEZIpIsIr8QkQwRSRORu7yeo4eIfCQi+SKyRUSeqFqjOc/+bgB2A6NEZDDwfeA2Vd2gquWquhu4AbhKRKZ7bfoRME1EQp3luUACkF6X160SQwpwGLjEWTXBiemLKut8gC1Vj6XTHPhbEVnnvA//EZFwr+Nzu4gcE5FsEXnY+7VFpJOIPCMiqc7tGRHp5Dz2hYjc4Ny/yHnNq5zlmSKyoz77Wdtnog77Mc3rc3jC+ezdC9wG/NL5/HzklD3XTFfLPp73s2VqZwnEnCMiMXh+dR90lqOBj4EngDDgAeBdEYkQkS54frFfoardgKlAvb5UVHUvcB9f/lIPqaFoLyAYiAbuBv7l9eX9L6DQKTPfudVlX0VELgJGAtuBmUCyqm6uEuMJYCNwudfqYuBD4GZn+Q7g5bq8bg1W82WyuARYA6ytsm6jqpbVsP2twF1AJBCA531CREYA/wvcDvQGegAxXts9DFyI5wfDWGAS8CvnsS+AGc796Xw1yU13Hm9qNe1HX+BT4B9AhBPvDlVdALwGPOV8fq6p5jnPt49w/s+WqYUlEAPwgYicBk4AGcCjzvpvA5+o6ieqWqmqy4F44Ern8Uo8v947q2qa84u9OZQBj6tqmap+AhQAQ8XTkXwD8KiqnlHVPcBLdXi+LCAHeA54UFVXAOFAWg3l0/B8cXl7GbhDRELwfKF+UM12Fzq/mM/eLqzh+b1rGxfjSSBrqqw73xf2C6p6QFWLgLf4sgZ5I7BEVVeragnwazzv2Vm34TmuGaqaCfwGT7I5G9PZWtclwB+8lmtLIDdV2e/e5ylbl/24FfhMVRc5n4FsVa3rj5Xz7SPU8Nmq43N3eJZADMB1Ti1iBjAMz5cpQF/gm95fBsA0IEpVC4Fv4alBpInIxyIyrJniy1bVcq/lM0BXPF/qfngS31ne92sSrqqhqjpcVf/urMsComooH+U8fo6qrnVe/2E8X9JF1Wy3UVVDvG4ba3j+1cAY55fvhXhqZPuAKGfdNM7f/+HddHb22ICn1nHueDjvWbZX2d7AMa/lY846gA3AEBHpieeL/GWgj9OsNKmWeN6qst+p5ylbl/3ow9f7l+rqfPsINX+2TB1YAjHnqOoXwIvAn51VJ4BXqnwZdFHVJ53yy1R1Np4v2H3A/znbFQJBXk/d63wv24iQM4Fyvtos06eBz/U5ni/ISd4rRaQPni/1VdVs8yrwCxrXfIWqHgZS8XTgH1fVAuehDc66rnia0eorDa/jISJBeJqxzkrF8yPhrFhnHap6BtgK/ARIVNVSYD3wc+CQqn4lodZBfT4TVZ0ABtbwWG2fnxr30TSeJRBT1TPAbBEZi+cL8hoRmSMiviIS6HQ8xohITxGZ5/SFlOCp+p9tHtkBXCKeaxaCOf9ZXSeBGPE6VbaunFNZ3wMeE5EgpwZ0Ry2b1fRcB/B06L8mnhMHfEVkJPAuni/Oz6rZ7O/AbJrm7Kg1eL6c13itW+usi6+hhlObd4CrnQ7oAOBxvvo/vwj4ldOnFQ48guc9P+sL4Id82Vy1qspyfdTnM1HVa8AsEblJRPzEc+LE2eatk8CA82xb2z6aRrAEYr7CaSd+GXjE6UCeB/wPnl/7J4D/wvO58cHz5ZaKpz9hOnC/8xzLgTfxnJm0FVhynpf8HM9ZR+kiUt9fteD5QgvG0/zxCp4vjJIGPM/Z53oOzxfMGSART5PHdapaWbWwquao6gptmkl1vsDTeex9BtkaZ12DEpTTJ/UD4HU8tZFTQLJXkSfw9GklALuAbc4675i6eb1+1eX6xFKfz0TVbY/j6Xf7BZ7P2g48HeIAzwMjnCbW6vqhattH0whiE0qZ9kRE/gj0UtU6nY1Vy3P9BrgeuERVcxsdnDHtjNVATJsmIsNEZIxzWu4kPKdivt8Uz62qjwIL8PSBGGOqsBqIadNEZCKeZqveeNrDFwBPNlGzkjHmPCyBGGOMaRBrwjLGGNMgfrUXaT/Cw8O1X79+bodhjDFtytatW7NUtepoDB0rgfTr14/4+Hi3wzDGmDZFRI5Vt96asIwxxjSIJRBjjDENYgnEGGNMg1gCMcYY0yCWQIwxxjSIJRBjjDENYgnEGGNMg1gCaSGf7krjUGZB7QWNMcZRWaks2nycnMJSt0OpliWQesg8XcI/ViRRWFJee2Ev6w5mcf9r23hscXNNGW6MaY9WJ2Xy0Hu7eODtndR33MLP9pxk2e702gs2giWQevj1B4n8ZfkB/rXyYJ23KSgp55fvJOAjsCYpi+PZZ5oxQmNMe7Jo83F8BD7fl8G721LqvF1KbhE/XLSNH72+vVlbPiyB1NHn+06ydHc6Ed068dyaI3VOBH/4ZC+peUX845YJ+Ai8seV4M0dqjGkPMvKL+WxvBndP68/EfqH85qPdpOcV12nbJz/dhyp08vfh1x8k1rv2UleWQOqgqLSCRz7czaDIrrx3/1T8fIXffbLnvNsUl1Xwt8+SeG3Tce6+qD9XjYnismE9eSs+mbKKr82OiqqycO0RLvvLKo5mFTbXrhhjWomH39/Fbc9tJON09Unh7a3JVFQqt0yK5akbx1JWUcm9r8SzP/30eZ9385EcPtqZyvemD+SXc4ex/lA2i3emNscuWAKpi398nkTyqSKeuG4UfcKC+MGlg1i2+ySrD2R+rWxeURkfbE9h7jOrefqzA1w1JooH5gwF4NbJfcgqKGHF3pNf2aasopL/eT+Rx5fs4XBmIb/+8Ku/GMqrSTjGmLbF+/941f4MXtt0nHUHs7nun+vYk5r/lbJnO8+nDOjBgIiu9A/vwtM3jeN4zhmu/PsaHv9oD/vS879WsyirqOTxJbuJCg7k/ukDuXVSLGNjgvntkj3kFZU1+T65OhqviMwF/gb4As+p6pNVHu8EvAxcAGQD31LVo85jD+GZvrQC+LGqLmuuOCsqlZviYrhwQA8A7p7Wn3e2JvP917ax4I4LmDownMSUPJ5atp/1B7Mor1QGhHfh5e9M4pIhX46APH1IJFHBgfzj84McyiykslLZd/I0246dIi2vmO/PGEh41048vmQPSxLSmD40gh8v2s7BjAL+87NLCAqo/e3al57PvrTTXDc+urkOhzEd3omcM6xJyuLWybF1Kv/vVQd5dtUh/nbLeKYM6MEjH+5mQEQX/vzNsXz/1W1c9+91jOsTwtiYYEK7BJCRX0LyqSJ+OXfYuee4YnQUkwf04E/
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADt0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjByYzMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy9h23ruAAAgAElEQVR4nOy9d5wcd33//3pvL7d7e10n3alasmS5YFs2NjYY3DAYMNWhBAgJIbSQQvmFEAIJIZhOiGkGvgEDMRhjSmyDK7Zxt2TZlmVJVtcVXW/b6+f3x2c+s5+Zndlyt3t70n2ej8c97m52duazszOfd39/iDEGhUKhUChkHM0egEKhUCiWHko4KBQKhaIEJRwUCoVCUYISDgqFQqEoQQkHhUKhUJSghINCoVAoSlDCQWELEd1PRO8t8/p3iejTizmm5Q4R/Z6I3l3H45X9jhXLFyUcThCI6AgRJYkoRkQjRPQjImpZxPP/BRE9JG9jjL2fMfa5BpxrLREx7bPGtM/+T9LrREQfJ6L92jU5RkT/SUQeaZ8face4xnTsr2vb/0L6XHnpXDEiut5mXPdr7z3LtP3X2vaX1/M6WMEYexVj7MfS2B+q9J5GUOk7Upz4KOFwYvFaxlgLgBcBOBvAJ5s8nkYT0T7v2wD8KxFdpW3/JoD3AXgXgBCAVwG4HMDPTe9/QdsHAEBELgDXAjho2u9RxliL9PPhMmMyH7MDwIUAxmv9cCcJ4jt6M4BPE9EVzR6Qoj4o4XACwhgbAXAnuJAAABDRBUT0CBHNENEzsharaZiHiChKRIeJ6B3a9s8S0U+l/YQ26JLPR0RbAHwXwIWaljijbf8REf2H9vfLiWiQiD5KRGNEdJyI3iMdo4OI/o+I5ojoSSL6j2q1XsbYowB2AzidiDYC+CCAdzDGHmWM5RhjuwG8CcDVRHSJ9Nb/A3AxEbVp/18F4FkAI9Wc14afAfgzInJq/78NwK8BZKTPej4RPap9F8eJ6HqTVXMlEe0jolki+jYRPSBcO8IaIKKvENG09n29Snrv/UT03jLficFNZLYuiOgKItqrnft6ACR/OCL6SyLao537TiJaU81FYYxtB/+O5HvS8lia5fd17T6ZI6JdRHS69tqPiLsr79bu1wfkMRDRS7T7Z1b7/RLTtfkcET2svfcuIurUXvMR0U+JaFL7Xp4koh7ttVYi+qH2XQ1p96b4fpctSjicgBBRH7i2fED7fxWA2wH8B4B2AB8D8Csi6iKiILim/SrGWAjASwA8Xcv5GGN7ALwfRQ07YrPrCgCtAFYB+CsA35Im5m8BiGv7vFv7qeazEhFdBGArgJ0ALgMwyBh7wjTGAQCPAbhS2pwC8FsAb9X+fxeAG6s5bxmGATwvncfqmHkA/wCgE9yquAxcoEGbrG4Bt/o6AOwD/05kXqxt7wTwJQA/JCLDJF7Dd6KjnftWAP+iHfsggIuk168B8M8A3gigC8CfANxU6bjaey8AcDqK92S5Y10J4GUANoHfL9cCmJQO9w4An9PG+DS4QAYRtYPf598Ev3ZfA3A7cetN8HYA7wHQDcAD/iwA/H5rBdCvvff9AJLaaz8CkANwCrhFfiWAZR+HUcLhxOI3RBQFMABgDMBntO1/DuAOxtgdjLECY+xuANsBvFp7vQCudfsZY8c1TbsRZAH8O2Msyxi7A0AMwKmaFvYmAJ9hjCUYY88D+HEVx5sAMAXgBwD+iTF2L/iEcdxm/+PgE5HMjQDeRUQRAJcA+I3F+y7QtEnxc0GFcYljbgZ3qzwqv8gY28EYe0yzao4A+J52boB/J7sZY7cyxnLgE53ZkjnKGPs+YywPfp16AfRUGFM1iHPfwhjLAviG6dzvB/AFxtgebWz/CeBFFayHCSJKAngUwLdRvL7ljpUFdwduBkDaPvJ3ejtj7EHGWBrAp8Cto34AVwPYzxj7iXZtbwKwF8Brpff+D2PsBcZYEsDNKFoyWXChcApjLK99R3Oa9fBqAH/PGIszxsYAfB1FhWLZooTDicXrNe3/5eAPVqe2fQ2At8gTHICLAfQyxuIA/gz8YT1ORLdrk1ojmNQmAkECQAv4hO0CF2oC+W87OhljbYyxLYyxb2rbJsAnSyt6tdd1GGMPaef/FIDbtEnDzGOMsYj081iFcd0K4FIAHwbwE/OLRLSJiG4jnjgwBz4xiu9qJaTPznjny0HTIUak1xPan/VIPrA6t/w9rAHwX9I9NAXudlpV5pid2tg+Cn5fuisdizF2H4Drwa3JMSK6gYjC0jHlMca0967Ufo6azn/UND5Z2In7D+Df050Afk5Ew0T0JSJya+N0gz8bYqzfA7c8ljVKOJyAMMYeADeFv6JtGgDwE9MEF2SMXaftfydj7ArwyXMvgO9r74sDCEiHXlHutAsY8ji42d4nbeuf57HuA9BPROfLGzXN8gIA91u856fgk9dCXUoA9An79wA+AAvhAOA74Nd5I2MsDO5eEW6h45Cug+Yu6is5QpVDsdhW7js9Dum6a+eWv4cBAH9juo/8jLFHyg6Ca+JfA3fjfbCaYzHGvskYOxfAaeDupY9Lh5TH2ALuKh3WfsxWzGoAQ+XGp50vyxj7N8bYaeBuvNeAuwQHAKTBFRExzjBjbGulY57sKOFw4vINAFcQT6v8KYDXEtEricipBd9eTkR9RNRDRNdosYc0uKunoB3jaQAvI6LVRNSK8tlPowD6SAqsVovmHrkVwGeJKKBZLu+q8Da7Y70AHoj9GfEgvJOItgL4FYBHANxj8bZvArgCwIPzOacN/wzgEs1tZCYEYA5ATPusH5Beux3AGUT0euKB/w+hvFAuh9V38jSAN2rX+RTw2I987q1E9Ebt3B8xnfu7AD6pXU8RqH1LDeO5DsAniMhX7lhEdB4RvVjT3OPgQqUgHefVRHSx9rk+B27ZDQC4A8AmIno7EbmI6M/AhcttlQZGRK8gojM0F+ccuJupoLmz7gLwVSIKE5GDiDaQMbFhWaKEwwkKY2wcXBP+V+3BEQHAcXBt6OPg368DwD+Ca11T4L7vD2jHuBvAL8AzeHag/EN2H3g2yggRTZTZz44PgwcER8C17ZvAhdV8+DB4HOKn4K6D58DdC69njBXMOzPGphhj92pulLrAGBvWXFZWfAw8MBoFt9J+Ib1vAsBbwAPNk+CT23bM71pYfSdfB8+cGgWPV/zM4tzXaefeCOBh6fVfA/giuOtlDvy66plSVXA7gGkAf13hWGHw6zIN/r1NAviydJz/BY+nTQE4FzymBsbYJLjG/1HtPZ8A8Brtc1ViBXgiwByAPQAeQNHqexd48Pp5bUy3wN51uWygOj4vCkXVENEXAaxgjC242peI/g3AGwC8jDE2s+DBLSJE5ACPObyDMfbHZo+n2RDRj8Cz0f6l2WNZ7ijLQbEoENFmIjpTS009H9zd8et6HJsx9hkAN4DHHJY8mvsvQkReFOMRlYLgCsWi4qq8i0JRF0LgrqSV4C6Pr4LXINQFxphly4slyoXgrhPhyni9TRaVQtE0lFtJoVAoFCUot5JCoVAoSmiaW0nLS78RvPKTAbiBMfZfWon8LwCsBXAEwLWMselyx+rs7GRr165t6HgVCoXiZGPHjh0TjDFzVwEATXQrEVEveAXvU0QUAk+lfD2AvwAwxRi7jngL4DbG2P9X7ljbtm1j27dvb/iYFQqF4mSCiHYwxrZZvdY0t5LW4+cp7e8oeO7xKvB8fdF358fgAkOhUCgUi8iSiDkQ0VrwboiPA+iRmnCNwKbhGBG9j4i2E9H28fHl2kpfoVAoGkPThYPWO+VX4F0R5+TXtIpWS78XY+wGxtg2xti2ri5Ll5lCoVAo5klThYPWW+VXAH7GGLtV2zyqxSNEXGKsWeNTKBSK5UrThIPWEfKHAPZoHR0Fv0NxIZh3o46FUgqFQqGojmZWSF8E4J0AdhGRWJnsn8Gbgt1MRH8F3pTr2iaNT6FQKJYtTRMOWkdLsnn5ssUci0KhUCiMND0g3UwOjMXwtbtfwEP7J6DaiCgUCkWRZd14b8/xOVx/334UGHD1mb34+rUvgse
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"err = optimfuncQMF(b)\n",
"print(err)\n",
"\n",
"#Restore symmetric upper half of window:\n",
"h = np.concatenate((xmin, np.flipud(xmin)))\n",
"plt.plot(h)\n",
"plt.title('Resulting PQMF Window Function')\n",
"plt.xlabel('Sample')\n",
"plt.ylabel('Value')\n",
"plt.show()\n",
"\n",
"f, H = sig.freqz(h)\n",
"plt.plot(f, 20*np.log10(np.abs(H)))\n",
"plt.title('Resulting PQMF Magnitude Response')\n",
"plt.xlabel('Normalized Frequency')\n",
"plt.ylabel('dB')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.7.4 64-bit ('base': conda)",
"language": "python",
"name": "python37464bitbaseconda58faf23c4b5f4fef93406f29a1005f35"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}