pySLFilters module
# This module contains all neccessary files to compute the filters used # in the pyShearLab2D toolbox. Most of these files are taken from different # MATLAB toolboxes and were translated to Python. # # Stefan Loock, February 2, 2017 [sloock@gwdg.de] import numpy as np from scipy import signal as signal ############################################################################## # # THIS IS A REWRITE OF THE ORIGINAL MATLAB IMPLEMENTATION OF # MakeONFilter.m FROM THE WaveLab850 TOOLBOX. -- Stefan Loock, Dec 2016. # # MakeONFilter -- Generate Orthonormal QMF Filter for Wavelet Transform # Usage # qmf = MakeONFilter(Type, Par) # Inputs # Type string, 'Haar', 'Beylkin', 'Coiflet', 'Daubechies', 'Symmlet', # 'Vaidyanathan', 'Battle' # # Outputs # qmf quadrature mirror filter # # Description # The Haar filter (which could be considered a Daubechies-2) was the # first wavelet, though not called as such, and is discontinuous. # # The Beylkin filter places roots for the frequency response function # close to the Nyquist frequency on the real axis. # # The Coiflet filters are designed to give both the mother and father # wavelets 2*Par vanishing moments; here Par may be one of 1,2,3,4 or 5. # # The Daubechies filters are minimal phase filters that generate wavelets # which have a minimal support for a given number of vanishing moments. # They are indexed by their length, Par, which may be one of # 4,6,8,10,12,14,16,18 or 20. The number of vanishing moments is par/2. # # Symmlets are also wavelets within a minimum size support for a given # number of vanishing moments, but they are as symmetrical as possible, # as opposed to the Daubechies filters which are highly asymmetrical. # They are indexed by Par, which specifies the number of vanishing # moments and is equal to half the size of the support. It ranges # from 4 to 10. # # The Vaidyanathan filter gives an exact reconstruction, but does not # satisfy any moment condition. The filter has been optimized for # speech coding. # # The Battle-Lemarie filter generate spline orthogonal wavelet basis. # The parameter Par gives the degree of the spline. The number of # vanishing moments is Par+1. # # See Also # FWT_PO, IWT_PO, FWT2_PO, IWT2_PO, WPAnalysis # # References # The books by Daubechies and Wickerhauser. # # Part of WaveLab850 (http://www-stat.stanford.edu/~wavelab/) def MakeONFilter(Type,Par=1): if Type == 'Haar': onFilter = np.array([1/np.sqrt(2), 1/np.sqrt(2)]) if Type == 'Beylkin': onFilter = np.array([.099305765374, .424215360813, .699825214057, .449718251149, -.110927598348, -.264497231446, .026900308804, .155538731877, -.017520746267, -.088543630623, .019679866044, .042916387274, -.017460408696, -.014365807969, .010040411845, .001484234782, -.002736031626, .000640485329]) if Type == 'Coiflet': if Par == 1: onFilter = np.array([.038580777748, -.126969125396, -.077161555496, .607491641386, .745687558934, .226584265197]) elif Par == 2: onFilter = np.array([.016387336463, -.041464936782, -.067372554722, .386110066823, .812723635450, .417005184424, -.076488599078, -.059434418646, .023680171947, .005611434819, -.001823208871, -.000720549445]) elif Par == 3: onFilter = np.array([-.003793512864, .007782596426, .023452696142, -.065771911281, -.061123390003, .405176902410, .793777222626, .428483476378, -.071799821619, -.082301927106, .034555027573, .015880544864, -.009007976137, -.002574517688, .001117518771, .000466216960, -.000070983303, -.000034599773]) elif Par == 4: onFilter = np.array([.000892313668, -.001629492013, -.007346166328, .016068943964, .026682300156, -.081266699680, -.056077313316, .415308407030, .782238930920, .434386056491, -.066627474263, -.096220442034, .039334427123, .025082261845, -.015211731527, -.005658286686, .003751436157, .001266561929, -.000589020757, -.000259974552, .000062339034, .000031229876, -.000003259680, -.000001784985]) elif Par == 5: onFilter = np.array([-.000212080863, .000358589677, .002178236305, -.004159358782, -.010131117538, .023408156762, .028168029062, -.091920010549, -.052043163216, .421566206729, .774289603740, .437991626228, -.062035963906, -.105574208706, .041289208741, .032683574283, -.019761779012, -.009164231153, .006764185419, .002433373209, -.001662863769, -.000638131296, .000302259520, .000140541149, -.000041340484, -.000021315014, .000003734597, .000002063806, -.000000167408, -.000000095158]) if Type == 'Daubechies': if Par == 4: onFilter = np.array([.482962913145, .836516303738, .224143868042, -.129409522551]) elif Par == 6: onFilter = np.array([.332670552950, .806891509311, .459877502118, -.135011020010, -.085441273882, .035226291882]) elif Par == 8: onFilter = np.array([.230377813309, .714846570553, .630880767930, -.027983769417, -.187034811719, .030841381836, .032883011667, -.010597401785]) elif Par == 10: onFilter = np.array([.160102397974, .603829269797, .724308528438, .138428145901, -.242294887066, -.032244869585, .077571493840, -.006241490213, -.012580751999, .003335725285]) elif Par == 12: onFilter = np.array([.111540743350, .494623890398, .751133908021, .315250351709, -.226264693965, -.129766867567, .097501605587, .027522865530, -.031582039317, .000553842201, .004777257511, -.001077301085]) elif Par == 14: onFilter = np.array([.077852054085, .396539319482, .729132090846, .469782287405, -.143906003929, -.224036184994, .071309219267, .080612609151, -.038029936935, -.016574541631, .012550998556, .000429577973, -.001801640704, .000353713800]) elif Par == 16: onFilter = np.array([.054415842243, .312871590914, .675630736297, .585354683654, -.015829105256, -.284015542962, .000472484574, .128747426620, -.017369301002, -.044088253931, .013981027917, .008746094047, -.004870352993, -.000391740373, .000675449406, -.000117476784]) elif Par==18: onFilter = np.array([.038077947364, .243834674613, .604823123690, .657288078051, .133197385825, -.293273783279, -.096840783223, .148540749338, .030725681479, -.067632829061, .000250947115, .022361662124, -.004723204758, -.004281503682, .001847646883, .000230385764, -.000251963189, .000039347320]) elif Par==20: onFilter = np.array([.026670057901, .188176800078, .527201188932, .688459039454, .281172343661, -.249846424327, -.195946274377, .127369340336, .093057364604, -.071394147166, -.029457536822, .033212674059, .003606553567, -.010733175483, .001395351747, .001992405295, -.000685856695, -.000116466855, .000093588670, -.000013264203]) if Type == 'Symmlet': if Par == 4: onFilter = np.array([-.107148901418, -.041910965125, .703739068656, 1.136658243408, .421234534204, -.140317624179, -.017824701442, .045570345896]) elif Par == 5: onFilter = np.array([.038654795955, .041746864422, -.055344186117, .281990696854, 1.023052966894, .896581648380, .023478923136, -.247951362613, -.029842499869, .027632152958]) elif Par == 6: onFilter = np.array([.021784700327, .004936612372, -.166863215412, -.068323121587, .694457972958, 1.113892783926, .477904371333, -.102724969862, -.029783751299, .063250562660, .002499922093, -.011031867509]) elif Par == 7: onFilter = np.array([.003792658534, -.001481225915, -.017870431651, .043155452582, .096014767936, -.070078291222, .024665659489, .758162601964, 1.085782709814, .408183939725, -.198056706807, -.152463871896, .005671342686, .014521394762]) elif Par == 8: onFilter = np.array([.002672793393, -.000428394300, -.021145686528, .005386388754, .069490465911, -.038493521263, -.073462508761, .515398670374, 1.099106630537, .680745347190, -.086653615406, -.202648655286, .010758611751, .044823623042, -.000766690896, -.004783458512]) elif Par == 9: onFilter = np.array([.001512487309, -.000669141509, -.014515578553, .012528896242, .087791251554, -.025786445930, -.270893783503, .049882830959, .873048407349, 1.015259790832, .337658923602, -.077172161097, .000825140929, .042744433602, -.016303351226, -.018769396836, .000876502539, .001981193736]) elif Par == 10: onFilter = np.array([.001089170447, .000135245020, -.012220642630, -.002072363923, .064950924579, .016418869426, -.225558972234, -.100240215031, .667071338154, 1.088251530500, .542813011213, -.050256540092, -.045240772218, .070703567550, .008152816799, -.028786231926, -.001137535314, .006495728375, .000080661204, -.000649589896]) if Type == 'Vaidyanathan': onFilter = np.array([-.000062906118, .000343631905, -.000453956620, -.000944897136, .002843834547, .000708137504, -.008839103409, .003153847056, .019687215010, -.014853448005, -.035470398607, .038742619293, .055892523691, -.077709750902, -.083928884366, .131971661417, .135084227129, -.194450471766, -.263494802488, .201612161775, .635601059872, .572797793211, .250184129505, .045799334111]) if Type == 'Battle': if Par == 1: onFilterTmp = np.array([0.578163, 0.280931, -0.0488618, -0.0367309, 0.012003, 0.00706442, -0.00274588, -0.00155701, 0.000652922, 0.000361781, -0.000158601, -0.0000867523]) elif Par == 3: onFilterTmp = np.array([0.541736, 0.30683, -0.035498, -0.0778079, 0.0226846, 0.0297468, -0.0121455, -0.0127154, 0.00614143, 0.00579932, -0.00307863, -0.00274529, 0.00154624, 0.00133086, -0.000780468, -0.00065562, 0.000395946, 0.000326749, -0.000201818, -0.000164264, 0.000103307]) elif Par == 5: onFilterTmp = np.array([0.528374, 0.312869, -0.0261771, -0.0914068, 0.0208414, 0.0433544, -0.0148537, -0.0229951, 0.00990635, 0.0128754, -0.00639886, -0.00746848, 0.00407882, 0.00444002, -0.00258816, -0.00268646, 0.00164132, 0.00164659, -0.00104207, -0.00101912, 0.000662836, 0.000635563, -0.000422485, -0.000398759, 0.000269842, 0.000251419, -0.000172685, -0.000159168, 0.000110709, 0.000101113]) onFilter = np.zeros(2*onFilterTmp.size-1) onFilter[onFilterTmp.size-1:2*onFilterTmp.size] = onFilterTmp; onFilter[0:onFilterTmp.size-1] = onFilterTmp[onFilterTmp.size-1:0:-1] return onFilter / np.linalg.norm(onFilter) # # Copyright (c) 1993-5. Jonathan Buckheit and David Donoho # # Part of Wavelab Version 850 # Built Tue Jan 3 13:20:40 EST 2006 # This is Copyrighted Material # For Copying permissions see COPYING.m # Comments? e-mail wavelab@stat.stanford.edu # ############################################################################## ############################################################################## # # THIS IS A REWRITE OF THE ORIGINAL MATLAB IMPLEMENTATION OF # dfilters.m FROM THE NONSUBSAMPLED CONTOURLET TOOLBOX. --Stefan Loock, Dec 2016. # THE FOLLOWING COMMENT IS FROM THE ORIGINAL FILE AND MAY ONLY APPLY IN PARTS # TO THIS IMPLEMENTATION: # DFILTERS Generate directional 2D filters # # [h0, h1] = dfilters(fname, type) # # Input: # fname: Filter name. Available 'fname' are: # 'haar': the "Haar" filters # 'vk': McClellan transformed of the filter from # the VK book # 'ko': orthogonal filter in the Kovacevic's paper # 'kos': smooth 'ko' filter # 'lax': 17 x 17 by Lu, Antoniou and Xu # 'sk': 9 x 9 by Shah and Kalker # 'cd': 7 and 9 McClellan transformed by # Cohen and Daubechies # 'pkva': ladder filters by Phong et al. # 'oqf_362': regular 3 x 6 filter # 'dvmlp' : regular linear phase biorthogonal filter with 3 dvm # 'sinc': ideal filter (*NO perfect recontruction*) # 'dmaxflat': diamond maxflat filters obtained from a three stage ladder # # type: 'd' or 'r' for decomposition or reconstruction filters # # Output: # h0, h1: diamond filter pair (lowpass and highpass) # To test those filters (for the PR condition for the FIR case), verify that: # conv2(h0, modulate2(h1, 'b')) + conv2(modulate2(h0, 'b'), h1) = 2 # (replace + with - for even size filters) # # To test for orthogonal filter # conv2(h, reverse2(h)) + modulate2(conv2(h, reverse2(h)), 'b') = 2 # # Part of the Nonsubsampled Contourlet Toolbox # (http://www.mathworks.de/matlabcentral/fileexchange/10049-nonsubsampled-contourlet-toolbox) def dfilters(fname, type): if fname == 'haar': if type.lower() == 'd': h0 = np.array([1, 1]) / np.sqrt(2) h1 = np.array([-1, 1]) / np.sqrt(2) else: h0 = np.array([1, 1]) / np.sqrt(2) h1 = np.array([1, -1]) / np.sqrt(2) elif fname == 'vk': # in Vetterli and Kovacevic book if type.lower() == 'd': h0 = np.array([1, 2, 1]) / 4 h1 = np.array([-1, -2, 6, -2, -1]) / 4 else: h0 = np.array([-1, 2, 6, 2, -1]) / 4 h1 = np.array([-1, 2, -1]) / 4 t = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]]) / 4 # diamon kernel h0 = mctrans(h0, t) h1 = mctrans(h1, t) elif fname == 'ko': # orthogonal filters in Kovacevics thesis a0 = 2 a1 = 0.5 a2 = 1 h0 = np.array([[0, -a1, -a0*a1, 0], [-a2, -a0*a2, -a0, 1], [0, a0*a1*a2, -a1*a2, 0]]) # h1 = qmf2(h0) h1 = np.array([[0, -a1*a2, -a0*a1*a2, 0], [1, a0, -a0*a2, a2], [0, -a0*a1, a1, 0]]) # normalize filter sum and norm norm = np.sqrt(2) / np.sum(h0) h0 = h0 * norm h1 = h1 * norm if type == 'r': # reverse filters for reconstruction h0 = h0[::-1] h1 = h1[::-1] elif fname == 'kos': # smooth orthogonal filters in Kovacevics thesis a0 = -np.sqrt(3) a1 = -np.sqrt(3) a2 = 2+np.sqrt(3) h0 = np.array([[0, -a1, -a0*a1, 0], [-a2, -a0*a2, -a0, 1], [0, a0*a1*a2, -a1*a2, 0]]) # h1 = qmf2(h0) h1 = np.array([[0, -a1*a2, -a0*a1*a2, 0], [1, a0, -a0*a2, a2], [0, -a0*a1, a1, 0]]) # normalize filter sum and norm norm = np.sqrt(2) / np.sum(h0) h0 = h0 * norm h1 = h1 * norm if type == 'r': # reverse filters for reconstruction h0 = h0[::-1] h1 = h1[::-1] elif fname == 'lax': # by lu, antoniou and xu h = np.array([[-1.2972901e-5, 1.2316237e-4, -7.5212207e-5, 6.3686104e-5, 9.4800610e-5, -7.5862919e-5, 2.9586164e-4, -1.8430337e-4], [1.2355540e-4, -1.2780882e-4, -1.9663685e-5, -4.5956538e-5, -6.5195193e-4, -2.4722942e-4, -2.1538331e-5, -7.0882131e-4], [-7.5319075e-5, -1.9350810e-5, -7.1947086e-4, 1.2295412e-3, 5.7411214e-4, 4.4705422e-4, 1.9623554e-3, 3.3596717e-4], [6.3400249e-5, -2.4947178e-4, 4.4905711e-4, -4.1053629e-3, -2.8588307e-3, 4.3782726e-3, -3.1690509e-3, -3.4371484e-3], [9.6404973e-5, -4.6116254e-5, 1.2371871e-3, -1.1675575e-2, 1.6173911e-2, -4.1197559e-3, 4.4911165e-3, 1.1635130e-2], [-7.6955555e-5, -6.5618379e-4, 5.7752252e-4, 1.6211426e-2, 2.1310378e-2, -2.8712621e-3, -4.8422645e-2, -5.9246338e-3], [2.9802986e-4, -2.1365364e-5, 1.9701350e-3, 4.5047673e-3, -4.8489158e-2, -3.1809526e-3, -2.9406153e-2, 1.8993868e-1], [-1.8556637e-4, -7.1279432e-4, 3.3839195e-4, 1.1662001e-2, -5.9398223e-3, -3.4467920e-3, 1.9006499e-1, 5.7235228e-1] ]) h0 = np.sqrt(2) * np.append(h, h[:,-2::-1], 1) h0 = np.append(h0, h0[-2::-1,:], 0) h1 = modulate2(h0, 'b') elif fname == 'sk': # by shah and kalker h = np.array([[0.621729, 0.161889, -0.0126949, -0.00542504, 0.00124838], [0.161889, -0.0353769, -0.0162751, -0.00499353, 0], [-0.0126949, -0.0162751, 0.00749029, 0, 0], [-0.00542504, 0.00499353, 0, 0, 0], [0.00124838, 0, 0, 0, 0]]) h0 = np.append(h[-1:0:-1, -1:0:-1], h[-1:0:-1,:], 1) h0 = np.append(h0, np.append(h[:,-1:0:-1], h, 1), 0)*np.sqrt(2) h1 = modulate2(h0, 'b') elif fname == 'dvmlp': q = np.sqrt(2) b = 0.02 b1 = b*b; h = np.array([[b/q, 0, -2*q*b, 0, 3*q*b, 0, -2*q*b, 0, b/q], [0, -1/(16*q), 0, 9/(16*q), 1/q, 9/(16*q), 0, -1/(16*q), 0], [b/q, 0, -2*q*b, 0, 3*q*b, 0, -2*q*b, 0, b/q]]) g0 = np.array([[-b1/q, 0, 4*b1*q, 0, -14*q*b1, 0, 28*q*b1, 0, -35*q*b1, 0, 28*q*b1, 0, -14*q*b1, 0, 4*b1*q, 0, -b1/q], [0, b/(8*q), 0, -13*b/(8*q), b/q, 33*b/(8*q), -2*q*b, -21*b/(8*q), 3*q*b, -21*b/(8*q), -2*q*b, 33*b/(8*q), b/q, -13*b/(8*q), 0, b/(8*q), 0], [-q*b1, 0, -1/(256*q) + 8*q*b1, 0, 9/(128*q) - 28*q*b1, -1/(q*16), -63/(256*q) + 56*q*b1, 9/(16*q), 87/(64*q)-70*q*b1, 9/(16*q), -63/(256*q) + 56*q*b1, -1/(q*16), 9/(128*q) - 28*q*b1, 0, -1/(256*q) + 8*q*b1, 0, -q*b1], [0, b/(8*q), 0, -13*b/(8*q), b/q, 33*b/(8*q), -2*q*b, -21*b/(8*q), 3*q*b, -21*b/(8*q), -2*q*b, 33*b/(8*q), b/q, -13*b/(8*q), 0, b/(8*q), 0], [-b1/q, 0, 4*b1*q, 0, -14*q*b1, 0, 28*q*b1, 0, -35*q*b1, 0, 28*q*b1, 0, -14*q*b1, 0, 4*b1*q, 0, -b1/q]]) h1 = modulate2(g0, 'b') h0 = h print(h1.shape) print(h0.shape) if type == 'r': h1 = modulate2(h, 'b') h0 = g0 elif fname == 'cd' or fname == '7-9': # by cohen and Daubechies h0 = np.array([0.026748757411, -0.016864118443, -0.078223266529, 0.266864118443, 0.602949018236, 0.266864118443, -0.078223266529, -0.016864118443, 0.026748757411]) g0 = np.array([-0.045635881557, -0.028771763114, 0.295635881557, 0.557543526229, 0.295635881557, -0.028771763114, -0.045635881557]) if type == 'd': h1 = modulate2(g0, 'c') else: h1 = modulate2(h0, 'c') h0 = g0 # use McClellan to obtain 2D filters t = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]])/4 # diamond kernel h0 = np.sqrt(2) * mctrans(h0, t) h1 = np.sqrt(2) * mctrans(h1, t) elif fname == 'oqf_362': h0 = np.sqrt(2) / 64 * np.array([[np.sqrt(15), -3, 0], [0, 5, np.sqrt(15)], [-2*np.sqrt(2), 30, 0], [0, 30, 2*np.sqrt(15)], [np.sqrt(15), 5, 0], [0, -3, -np.sqrt(15)]]) h1 = -modulate2(h0, 'b') h1 = -h1[::-1] if type == 'r': h0 = h0[::-1] h1 = -modulate2(h0, 'b') h1 = -h1[::-1] elif fname == 'test': h0 = np.array([[0, 1, 0], [1, 4, 1], [0, 1, 0]]) h1 = np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]]) elif fname == 'testDVM': h0 = np.array([[1, 1], [1, 1]]) / np.sqrt(2) h1 = np.array([[-1, 1], [1, -1]]) / np.sqrt(2) elif fname == 'qmf': # by Lu, antoniou and xu # ideal response window m = 2 n = 2 w1d = np.kaiser(4*m+1, 2.6) w = np.zeros((n+m+1,n+m+1)) for n1 in np.arange(-m,m+1): for n2 in np.arange(-n,n+1): w[n1+m,n2+n] = w1d[2*m+n1+n2]*w1d[2*m+n1-n2] h = np.zeros((n+m+1,n+m+1)) for n1 in np.arange(-m,m+1): for n2 in np.arange(-n,n+1): h[n1+m, n2+n] = 0.5*np.sinc((n1+n2)/2) * 0.5*np.sinc((n1-n2)/2) c = np.sum(h) h = np.sqrt(2) * h h0 = h * w h1 = modulate2(h0, 'b') elif fname == 'qmf2': # by Lu, Antoniou and Xu # ideal response window h = np.array([ [-0.001104, 0.002494, -0.001744, 0.004895, -0.000048, -0.000311], [0.008918, -0.002844, -0.025197, -0.017135, 0.003905, -0.000081], [-0.007587, -0.065904, 00.100431, -0.055878, 0.007023, 0.001504], [0.001725, 0.184162, 0.632115, 0.099414, -0.027006, -0.001110], [-0.017935, -0.000491, 0.191397, -0.001787, -0.010587, 0.002060], [0.001353, 0.005635, -0.001231, -0.009052, -0.002668, 0.000596]]) h0 = h/np.sum(h) h1 = modulate2(h0, 'b') elif fname == 'dmaxflat4': M1 = 1/np.sqrt(2) M2 = np.copy(M1) k1 = 1-np.sqrt(2) k3 = np.copy(k1) k2 = np.copy(M1) h = np.array([0.25*k2*k3, 0.5*k2, 1+0.5*k2*k3])*M1 h = np.append(h, h[-2::-1]) g = np.array([-0.125*k1*k2*k3, 0.25*k1*k2, -0.5*k1-0.5*k3-0.375*k1*k2*k3, 1+0.5*k1*k2])*M2 g = np.append(g, h[-2::-1]) B = dmaxflat(4,0) h0 = mctrans(h,B) g0 = mctrans(g,B) h0 = np.sqrt(2) * h0 / np.sum(h0) g0 = np.sqrt(2) * g0 / np.sum(g0) h1 = modulate2(g0, 'b') if type == 'r': h1 = modulate2(h0, 'b') h0 = g0 elif fname == 'dmaxflat5': M1 = 1/np.sqrt(2) M2 = M1 k1 = 1-np.sqrt(2) k3 = k1 k2 = M1 h = np.array([0.25*k2*k3, 0.5*k2, 1+0.5*k2*k3])*M1 h = np.append(h, h[-2::-1]) g = np.array([-0.125*k1*k2*k3, 0.25*k1*k2, -0.5*k1-0.5*k3-0.375*k1*k2*k3, 1+0.5*k1*k2])*M2 g = np.append(g, h[-2::-1]) B = dmaxflat(5,0) h0 = mctrans(h,B) g0 = mctrans(g,B) h0 = np.sqrt(2) * h0 / np.sum(h0) g0 = np.sqrt(2) * g0 / np.sum(g0) h1 = modulate2(g0, 'b') if type == 'r': h1 = modulate2(h0, 'b') h0 = g0 elif fname == 'dmaxflat6': M1 = 1/np.sqrt(2) M2 = M1 k1 = 1-np.sqrt(2) k3 = k1 k2 = M1 h = np.array([0.25*k2*k3, 0.5*k2, 1+0.5*k2*k3])*M1 h = np.append(h, h[-2::-1]) g = np.array([-0.125*k1*k2*k3, 0.25*k1*k2, -0.5*k1-0.5*k3-0.375*k1*k2*k3, 1+0.5*k1*k2])*M2 g = np.append(g, h[-2::-1]) B = dmaxflat(6,0) h0 = mctrans(h,B) g0 = mctrans(g,B) h0 = np.sqrt(2) * h0 / np.sum(h0) g0 = np.sqrt(2) * g0 / np.sum(g0) h1 = modulate2(g0, 'b') if type == 'r': h1 = modulate2(h0, 'b') h0 = g0 elif fname == 'dmaxflat7': M1 = 1/np.sqrt(2) M2 = M1 k1 = 1-np.sqrt(2) k3 = k1 k2 = M1 h = np.array([0.25*k2*k3, 0.5*k2, 1+0.5*k2*k3])*M1 h = np.append(h, h[-2::-1]) g = np.array([-0.125*k1*k2*k3, 0.25*k1*k2, -0.5*k1-0.5*k3-0.375*k1*k2*k3, 1+0.5*k1*k2])*M2 g = np.append(g, h[-2::-1]) B = dmaxflat(7,0) h0 = mctrans(h,B) g0 = mctrans(g,B) h0 = np.sqrt(2) * h0 / np.sum(h0) g0 = np.sqrt(2) * g0 / np.sum(g0) h1 = modulate2(g0, 'b') if type == 'r': h1 = modulate2(h0, 'b') h0 = g0 # The original file supports a case "otherwise" for unrecognized filters # and computes simple 1D wavelet filters for them using wfilters.m # I think we don't need this and skip this for the time being. # IN ORIGINAL MATLAB VERSION: # otherwise # % Assume the "degenerated" case: 1D wavelet filters # [h0,h1] = wfilters(fname, type); return h0, h1 # ############################################################################## ############################################################################## # # THIS IS A REWRITE OF THE ORIGINAL MATLAB IMPLEMENTATION OF dmaxflat.m # FROM THE Nonsubsampled Contourlet Toolbox. -- Stefan Loock, Dec 2016. # # returns 2-D diamond maxflat filters of order 'N' # the filters are nonseparable and 'd' is the (0,0) coefficient, being 1 or 0 # depending on use. # by Arthur L. da Cunha, University of Illinois Urbana-Champaign # Aug 2004 def dmaxflat(N,d): if (N > 7) or (N < 1): print('Error: N must be in {1,2,...,7}') return 0 if N == 1: h = np.array([[0, 1, 0],[1, 0, 1],[0, 1, 0]])/4 h[1,1] = d elif N == 2: h = np.array([[0, -1, 0],[-1, 0, 10], [0, 10, 0]]) h = np.append(h, np.fliplr(h[:,0:-1]), 1) h = np.append(h, np.flipud(h[0:-1,:]), 0)/32 h[2,2] = d elif N == 3: h = np.array([[0, 3, 0, 2],[3, 0, -27, 0],[0, -27, 0, 174], [2, 0, 174, 0]]) h = np.append(h, np.fliplr(h[:, 0:-1]), 1) h = np.append(h, np.flipud(h[0:-1,:]),0) h[3,3] = d elif N == 4: h = np.array([[0, -5, 0, -3, 0], [-5, 0, 52, 0, 34], [0, 52, 0, -276, 0], [-3, 0, -276, 0, 1454], [0, 34, 0, 1454, 0]])/np.power(2,12) h = np.append(h, np.fliplr(h[:,0:-1]),1) h = np.append(h, np.flipud(h[0:-1,:]),0) h[4,4] = d elif N == 5: h = np.array([[0, 35, 0, 20, 0, 18], [35, 0, -425, 0, -250, 0], [0, -425, 0, 2500, 0, 1610], [20, 0, 2500, 0, -10200, 0], [0, -250, 0, -10200, 0, 47780], [18, 0, 1610, 0, 47780, 0]])/np.power(2,17) h = np.append(h, np.fliplr(h[:,0:-1]),1) h = np.append(h, np.flipud(h[0:-1,:]),0) h[5,5] = d elif N == 6: h = np.array([[0, -63, 0, -35, 0, -30, 0], [-63, 0, 882, 0, 495, 0, 444], [0, 882, 0, -5910, 0, -3420, 0], [-35, 0, -5910, 0, 25875, 0, 16460], [0, 495, 0, 25875, 0, -89730, 0], [-30, 0, -3420, 0, -89730, 0, 389112], [0, 44, 0, 16460, 0, 389112, 0]])/np.power(2,20) h = np.append(h, np.fliplr(h[:,0:-1]),1) h = np.append(h, np.flipud(h[0:-1,:]),0) h[6,6] = d elif N == 7: h = np.array([[0, 231, 0, 126, 0, 105, 0, 100], [231, 0, -3675, 0, -2009, 0, -1715, 0], [0, -3675, 0, 27930, 0, 15435, 0, 13804], [126, 0, 27930, 0, -136514, 0, -77910, 0], [0, -2009, 0, -136514, 0, 495145, 0, 311780], [105, 0, 15435, 0, 495145, 0, -1535709, 0], [0, -1715, 0, -77910, 0, -1534709, 0, 6305740], [100, 0, 13804, 0, 311780, 0, 6305740, 0]])/np.power(2,24) h = np.append(h, np.fliplr(h[:,0:-1]),1) h = np.append(h, np.flipud(h[0:-1,:]),0) h[7,7] = d return h # ############################################################################## ############################################################################## # # THIS IS A REWRITE OF THE ORIGINAL MATLAB IMPLEMENTATION OF mctrans.m # FROM THE Nonsubsampled Contourlet Toolbox. -- Stefan Loock, Dec 2016. # MCTRANS McClellan transformation # # # H = mctrans(B,T) produces the 2-D FIR filter H that # corresponds to the 1-D FIR filter B using the transform T. # # # Convert the 1-D filter b to SUM_n a(n) cos(wn) form # # Part of the Nonsubsampled Contourlet Toolbox # (http://www.mathworks.de/matlabcentral/fileexchange/10049-nonsubsampled-contourlet-toolbox) def mctrans(b,t): # Convert the 1-D filter b to SUM_n a(n) cos(wn) form # if mod(n,2) != 0 -> error n = (b.size-1)//2 b = np.fft.fftshift(b[::-1]) #inverse fftshift b = b[::-1] a = np.zeros(n+1) a[0] = b[0] a[1:n+1] = 2*b[1:n+1] inset = np.floor((np.asarray(t.shape)-1)/2) inset = inset.astype(int) # Use Chebyshev polynomials to compute h P0 = 1 P1 = t; h = a[1]*P1; rows = int(inset[0]+1) cols = int(inset[1]+1) h[rows-1,cols-1] = h[rows-1,cols-1]+a[0]*P0; for i in range(3,n+2): P2 = 2*signal.convolve2d(t, P1) rows = (rows + inset[0]).astype(int) cols = (cols + inset[1]).astype(int) if i == 3: P2[rows-1,cols-1] = P2[rows-1,cols-1] - P0 else: P2[rows[0]-1:rows[-1],cols[0]-1:cols[-1]] = P2[rows[0]-1:rows[-1], cols[0]-1:cols[-1]] - P0 rows = inset[0] + np.arange(np.asarray(P1.shape)[0])+1 rows = rows.astype(int) cols = inset[1] + np.arange(np.asarray(P1.shape)[1])+1 cols = cols.astype(int) hh = h; h = a[i-1]*P2 h[rows[0]-1:rows[-1], cols[0]-1:cols[-1]] = h[rows[0]-1:rows[-1], cols[0]-1:cols[-1]] + hh P0 = P1; P1 = P2; h = np.rot90(h,2) return h # ############################################################################## ############################################################################## # # THIS IS A REWRITE OF THE ORIGINAL MATLAB IMPLEMENTATION OF # modulate2.m FROM THE Nonsubsampled Contourlet Toolbox. # -- Stefan Loock, Dec 2016. # # MODULATE2 2D modulation # # y = modulate2(x, type, [center]) # # With TYPE = {'r', 'c' or 'b'} for modulate along the row, or column or # both directions. # # CENTER secify the origin of modulation as floor(size(x)/2)+1+center # (default is [0, 0]) # # Part of the Nonsubsampled Contourlet Toolbox # (http://www.mathworks.de/matlabcentral/fileexchange/10049-nonsubsampled-contourlet-toolbox) def modulate2(x, type, center=np.array([0, 0])): size = np.asarray(x.shape) if x.ndim == 1: if np.array_equal(center, [0, 0]): center = 0 origin = np.floor(size/2)+1+center n1 = np.arange(size[0])-origin[0]+1 if x.ndim == 2: n2 = np.arange(size[1])-origin[1]+1 else: n2 = n1 if type == 'r': m1 = np.power(-1,n1) if x.ndim == 1: y = x*m1 else: y = x * np.transpose(np.tile(m1, (size[1], 1))) elif type == 'c': m2 = np.power(-1,n2) if x.ndim == 1: y = x*m2 else: y = x * np.tile(m2, np.array([size[0], 1])) elif type == 'b': m1 = np.power(-1,n1) m2 = np.power(-1,n2) m = np.outer(m1, m2) if x.ndim == 1: y = x * m1 else: y = x * m return y # ############################################################################## ############################################################################## # # THIS IS A REWRITE OF THE ORIGINAL MATLAB IMPLEMENTATION OF # MirrorFilt.m FROM THE WaveLab850 TOOLBOX. -- Stefan Loock, Dec 2016. # # MirrorFilt -- Apply (-1)^t modulation # Usage # h = MirrorFilt(l) # Inputs # l 1-d signal # Outputs # h 1-d signal with DC frequency content shifted # to Nyquist frequency # # Description # h(t) = (-1)^(t-1) * x(t), 1 <= t <= length(x) # # See Also # DyadDownHi # # Part of WaveLab850 (http://www-stat.stanford.edu/~wavelab/) import numpy as np def MirrorFilt(x): return np.power(-1,np.arange(x.size))*x; # # Copyright (c) 1993. Iain M. Johnstone # # Part of Wavelab Version 850 # Built Tue Jan 3 13:20:40 EST 2006 # This is Copyrighted Material # For Copying permissions see COPYING.m # Comments? e-mail wavelab@stat.stanford.edu # ##############################################################################
Functions
def MakeONFilter(
Type, Par=1)
def MakeONFilter(Type,Par=1): if Type == 'Haar': onFilter = np.array([1/np.sqrt(2), 1/np.sqrt(2)]) if Type == 'Beylkin': onFilter = np.array([.099305765374, .424215360813, .699825214057, .449718251149, -.110927598348, -.264497231446, .026900308804, .155538731877, -.017520746267, -.088543630623, .019679866044, .042916387274, -.017460408696, -.014365807969, .010040411845, .001484234782, -.002736031626, .000640485329]) if Type == 'Coiflet': if Par == 1: onFilter = np.array([.038580777748, -.126969125396, -.077161555496, .607491641386, .745687558934, .226584265197]) elif Par == 2: onFilter = np.array([.016387336463, -.041464936782, -.067372554722, .386110066823, .812723635450, .417005184424, -.076488599078, -.059434418646, .023680171947, .005611434819, -.001823208871, -.000720549445]) elif Par == 3: onFilter = np.array([-.003793512864, .007782596426, .023452696142, -.065771911281, -.061123390003, .405176902410, .793777222626, .428483476378, -.071799821619, -.082301927106, .034555027573, .015880544864, -.009007976137, -.002574517688, .001117518771, .000466216960, -.000070983303, -.000034599773]) elif Par == 4: onFilter = np.array([.000892313668, -.001629492013, -.007346166328, .016068943964, .026682300156, -.081266699680, -.056077313316, .415308407030, .782238930920, .434386056491, -.066627474263, -.096220442034, .039334427123, .025082261845, -.015211731527, -.005658286686, .003751436157, .001266561929, -.000589020757, -.000259974552, .000062339034, .000031229876, -.000003259680, -.000001784985]) elif Par == 5: onFilter = np.array([-.000212080863, .000358589677, .002178236305, -.004159358782, -.010131117538, .023408156762, .028168029062, -.091920010549, -.052043163216, .421566206729, .774289603740, .437991626228, -.062035963906, -.105574208706, .041289208741, .032683574283, -.019761779012, -.009164231153, .006764185419, .002433373209, -.001662863769, -.000638131296, .000302259520, .000140541149, -.000041340484, -.000021315014, .000003734597, .000002063806, -.000000167408, -.000000095158]) if Type == 'Daubechies': if Par == 4: onFilter = np.array([.482962913145, .836516303738, .224143868042, -.129409522551]) elif Par == 6: onFilter = np.array([.332670552950, .806891509311, .459877502118, -.135011020010, -.085441273882, .035226291882]) elif Par == 8: onFilter = np.array([.230377813309, .714846570553, .630880767930, -.027983769417, -.187034811719, .030841381836, .032883011667, -.010597401785]) elif Par == 10: onFilter = np.array([.160102397974, .603829269797, .724308528438, .138428145901, -.242294887066, -.032244869585, .077571493840, -.006241490213, -.012580751999, .003335725285]) elif Par == 12: onFilter = np.array([.111540743350, .494623890398, .751133908021, .315250351709, -.226264693965, -.129766867567, .097501605587, .027522865530, -.031582039317, .000553842201, .004777257511, -.001077301085]) elif Par == 14: onFilter = np.array([.077852054085, .396539319482, .729132090846, .469782287405, -.143906003929, -.224036184994, .071309219267, .080612609151, -.038029936935, -.016574541631, .012550998556, .000429577973, -.001801640704, .000353713800]) elif Par == 16: onFilter = np.array([.054415842243, .312871590914, .675630736297, .585354683654, -.015829105256, -.284015542962, .000472484574, .128747426620, -.017369301002, -.044088253931, .013981027917, .008746094047, -.004870352993, -.000391740373, .000675449406, -.000117476784]) elif Par==18: onFilter = np.array([.038077947364, .243834674613, .604823123690, .657288078051, .133197385825, -.293273783279, -.096840783223, .148540749338, .030725681479, -.067632829061, .000250947115, .022361662124, -.004723204758, -.004281503682, .001847646883, .000230385764, -.000251963189, .000039347320]) elif Par==20: onFilter = np.array([.026670057901, .188176800078, .527201188932, .688459039454, .281172343661, -.249846424327, -.195946274377, .127369340336, .093057364604, -.071394147166, -.029457536822, .033212674059, .003606553567, -.010733175483, .001395351747, .001992405295, -.000685856695, -.000116466855, .000093588670, -.000013264203]) if Type == 'Symmlet': if Par == 4: onFilter = np.array([-.107148901418, -.041910965125, .703739068656, 1.136658243408, .421234534204, -.140317624179, -.017824701442, .045570345896]) elif Par == 5: onFilter = np.array([.038654795955, .041746864422, -.055344186117, .281990696854, 1.023052966894, .896581648380, .023478923136, -.247951362613, -.029842499869, .027632152958]) elif Par == 6: onFilter = np.array([.021784700327, .004936612372, -.166863215412, -.068323121587, .694457972958, 1.113892783926, .477904371333, -.102724969862, -.029783751299, .063250562660, .002499922093, -.011031867509]) elif Par == 7: onFilter = np.array([.003792658534, -.001481225915, -.017870431651, .043155452582, .096014767936, -.070078291222, .024665659489, .758162601964, 1.085782709814, .408183939725, -.198056706807, -.152463871896, .005671342686, .014521394762]) elif Par == 8: onFilter = np.array([.002672793393, -.000428394300, -.021145686528, .005386388754, .069490465911, -.038493521263, -.073462508761, .515398670374, 1.099106630537, .680745347190, -.086653615406, -.202648655286, .010758611751, .044823623042, -.000766690896, -.004783458512]) elif Par == 9: onFilter = np.array([.001512487309, -.000669141509, -.014515578553, .012528896242, .087791251554, -.025786445930, -.270893783503, .049882830959, .873048407349, 1.015259790832, .337658923602, -.077172161097, .000825140929, .042744433602, -.016303351226, -.018769396836, .000876502539, .001981193736]) elif Par == 10: onFilter = np.array([.001089170447, .000135245020, -.012220642630, -.002072363923, .064950924579, .016418869426, -.225558972234, -.100240215031, .667071338154, 1.088251530500, .542813011213, -.050256540092, -.045240772218, .070703567550, .008152816799, -.028786231926, -.001137535314, .006495728375, .000080661204, -.000649589896]) if Type == 'Vaidyanathan': onFilter = np.array([-.000062906118, .000343631905, -.000453956620, -.000944897136, .002843834547, .000708137504, -.008839103409, .003153847056, .019687215010, -.014853448005, -.035470398607, .038742619293, .055892523691, -.077709750902, -.083928884366, .131971661417, .135084227129, -.194450471766, -.263494802488, .201612161775, .635601059872, .572797793211, .250184129505, .045799334111]) if Type == 'Battle': if Par == 1: onFilterTmp = np.array([0.578163, 0.280931, -0.0488618, -0.0367309, 0.012003, 0.00706442, -0.00274588, -0.00155701, 0.000652922, 0.000361781, -0.000158601, -0.0000867523]) elif Par == 3: onFilterTmp = np.array([0.541736, 0.30683, -0.035498, -0.0778079, 0.0226846, 0.0297468, -0.0121455, -0.0127154, 0.00614143, 0.00579932, -0.00307863, -0.00274529, 0.00154624, 0.00133086, -0.000780468, -0.00065562, 0.000395946, 0.000326749, -0.000201818, -0.000164264, 0.000103307]) elif Par == 5: onFilterTmp = np.array([0.528374, 0.312869, -0.0261771, -0.0914068, 0.0208414, 0.0433544, -0.0148537, -0.0229951, 0.00990635, 0.0128754, -0.00639886, -0.00746848, 0.00407882, 0.00444002, -0.00258816, -0.00268646, 0.00164132, 0.00164659, -0.00104207, -0.00101912, 0.000662836, 0.000635563, -0.000422485, -0.000398759, 0.000269842, 0.000251419, -0.000172685, -0.000159168, 0.000110709, 0.000101113]) onFilter = np.zeros(2*onFilterTmp.size-1) onFilter[onFilterTmp.size-1:2*onFilterTmp.size] = onFilterTmp; onFilter[0:onFilterTmp.size-1] = onFilterTmp[onFilterTmp.size-1:0:-1] return onFilter / np.linalg.norm(onFilter)
def MirrorFilt(
x)
def MirrorFilt(x): return np.power(-1,np.arange(x.size))*x;
def dfilters(
fname, type)
def dfilters(fname, type): if fname == 'haar': if type.lower() == 'd': h0 = np.array([1, 1]) / np.sqrt(2) h1 = np.array([-1, 1]) / np.sqrt(2) else: h0 = np.array([1, 1]) / np.sqrt(2) h1 = np.array([1, -1]) / np.sqrt(2) elif fname == 'vk': # in Vetterli and Kovacevic book if type.lower() == 'd': h0 = np.array([1, 2, 1]) / 4 h1 = np.array([-1, -2, 6, -2, -1]) / 4 else: h0 = np.array([-1, 2, 6, 2, -1]) / 4 h1 = np.array([-1, 2, -1]) / 4 t = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]]) / 4 # diamon kernel h0 = mctrans(h0, t) h1 = mctrans(h1, t) elif fname == 'ko': # orthogonal filters in Kovacevics thesis a0 = 2 a1 = 0.5 a2 = 1 h0 = np.array([[0, -a1, -a0*a1, 0], [-a2, -a0*a2, -a0, 1], [0, a0*a1*a2, -a1*a2, 0]]) # h1 = qmf2(h0) h1 = np.array([[0, -a1*a2, -a0*a1*a2, 0], [1, a0, -a0*a2, a2], [0, -a0*a1, a1, 0]]) # normalize filter sum and norm norm = np.sqrt(2) / np.sum(h0) h0 = h0 * norm h1 = h1 * norm if type == 'r': # reverse filters for reconstruction h0 = h0[::-1] h1 = h1[::-1] elif fname == 'kos': # smooth orthogonal filters in Kovacevics thesis a0 = -np.sqrt(3) a1 = -np.sqrt(3) a2 = 2+np.sqrt(3) h0 = np.array([[0, -a1, -a0*a1, 0], [-a2, -a0*a2, -a0, 1], [0, a0*a1*a2, -a1*a2, 0]]) # h1 = qmf2(h0) h1 = np.array([[0, -a1*a2, -a0*a1*a2, 0], [1, a0, -a0*a2, a2], [0, -a0*a1, a1, 0]]) # normalize filter sum and norm norm = np.sqrt(2) / np.sum(h0) h0 = h0 * norm h1 = h1 * norm if type == 'r': # reverse filters for reconstruction h0 = h0[::-1] h1 = h1[::-1] elif fname == 'lax': # by lu, antoniou and xu h = np.array([[-1.2972901e-5, 1.2316237e-4, -7.5212207e-5, 6.3686104e-5, 9.4800610e-5, -7.5862919e-5, 2.9586164e-4, -1.8430337e-4], [1.2355540e-4, -1.2780882e-4, -1.9663685e-5, -4.5956538e-5, -6.5195193e-4, -2.4722942e-4, -2.1538331e-5, -7.0882131e-4], [-7.5319075e-5, -1.9350810e-5, -7.1947086e-4, 1.2295412e-3, 5.7411214e-4, 4.4705422e-4, 1.9623554e-3, 3.3596717e-4], [6.3400249e-5, -2.4947178e-4, 4.4905711e-4, -4.1053629e-3, -2.8588307e-3, 4.3782726e-3, -3.1690509e-3, -3.4371484e-3], [9.6404973e-5, -4.6116254e-5, 1.2371871e-3, -1.1675575e-2, 1.6173911e-2, -4.1197559e-3, 4.4911165e-3, 1.1635130e-2], [-7.6955555e-5, -6.5618379e-4, 5.7752252e-4, 1.6211426e-2, 2.1310378e-2, -2.8712621e-3, -4.8422645e-2, -5.9246338e-3], [2.9802986e-4, -2.1365364e-5, 1.9701350e-3, 4.5047673e-3, -4.8489158e-2, -3.1809526e-3, -2.9406153e-2, 1.8993868e-1], [-1.8556637e-4, -7.1279432e-4, 3.3839195e-4, 1.1662001e-2, -5.9398223e-3, -3.4467920e-3, 1.9006499e-1, 5.7235228e-1] ]) h0 = np.sqrt(2) * np.append(h, h[:,-2::-1], 1) h0 = np.append(h0, h0[-2::-1,:], 0) h1 = modulate2(h0, 'b') elif fname == 'sk': # by shah and kalker h = np.array([[0.621729, 0.161889, -0.0126949, -0.00542504, 0.00124838], [0.161889, -0.0353769, -0.0162751, -0.00499353, 0], [-0.0126949, -0.0162751, 0.00749029, 0, 0], [-0.00542504, 0.00499353, 0, 0, 0], [0.00124838, 0, 0, 0, 0]]) h0 = np.append(h[-1:0:-1, -1:0:-1], h[-1:0:-1,:], 1) h0 = np.append(h0, np.append(h[:,-1:0:-1], h, 1), 0)*np.sqrt(2) h1 = modulate2(h0, 'b') elif fname == 'dvmlp': q = np.sqrt(2) b = 0.02 b1 = b*b; h = np.array([[b/q, 0, -2*q*b, 0, 3*q*b, 0, -2*q*b, 0, b/q], [0, -1/(16*q), 0, 9/(16*q), 1/q, 9/(16*q), 0, -1/(16*q), 0], [b/q, 0, -2*q*b, 0, 3*q*b, 0, -2*q*b, 0, b/q]]) g0 = np.array([[-b1/q, 0, 4*b1*q, 0, -14*q*b1, 0, 28*q*b1, 0, -35*q*b1, 0, 28*q*b1, 0, -14*q*b1, 0, 4*b1*q, 0, -b1/q], [0, b/(8*q), 0, -13*b/(8*q), b/q, 33*b/(8*q), -2*q*b, -21*b/(8*q), 3*q*b, -21*b/(8*q), -2*q*b, 33*b/(8*q), b/q, -13*b/(8*q), 0, b/(8*q), 0], [-q*b1, 0, -1/(256*q) + 8*q*b1, 0, 9/(128*q) - 28*q*b1, -1/(q*16), -63/(256*q) + 56*q*b1, 9/(16*q), 87/(64*q)-70*q*b1, 9/(16*q), -63/(256*q) + 56*q*b1, -1/(q*16), 9/(128*q) - 28*q*b1, 0, -1/(256*q) + 8*q*b1, 0, -q*b1], [0, b/(8*q), 0, -13*b/(8*q), b/q, 33*b/(8*q), -2*q*b, -21*b/(8*q), 3*q*b, -21*b/(8*q), -2*q*b, 33*b/(8*q), b/q, -13*b/(8*q), 0, b/(8*q), 0], [-b1/q, 0, 4*b1*q, 0, -14*q*b1, 0, 28*q*b1, 0, -35*q*b1, 0, 28*q*b1, 0, -14*q*b1, 0, 4*b1*q, 0, -b1/q]]) h1 = modulate2(g0, 'b') h0 = h print(h1.shape) print(h0.shape) if type == 'r': h1 = modulate2(h, 'b') h0 = g0 elif fname == 'cd' or fname == '7-9': # by cohen and Daubechies h0 = np.array([0.026748757411, -0.016864118443, -0.078223266529, 0.266864118443, 0.602949018236, 0.266864118443, -0.078223266529, -0.016864118443, 0.026748757411]) g0 = np.array([-0.045635881557, -0.028771763114, 0.295635881557, 0.557543526229, 0.295635881557, -0.028771763114, -0.045635881557]) if type == 'd': h1 = modulate2(g0, 'c') else: h1 = modulate2(h0, 'c') h0 = g0 # use McClellan to obtain 2D filters t = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]])/4 # diamond kernel h0 = np.sqrt(2) * mctrans(h0, t) h1 = np.sqrt(2) * mctrans(h1, t) elif fname == 'oqf_362': h0 = np.sqrt(2) / 64 * np.array([[np.sqrt(15), -3, 0], [0, 5, np.sqrt(15)], [-2*np.sqrt(2), 30, 0], [0, 30, 2*np.sqrt(15)], [np.sqrt(15), 5, 0], [0, -3, -np.sqrt(15)]]) h1 = -modulate2(h0, 'b') h1 = -h1[::-1] if type == 'r': h0 = h0[::-1] h1 = -modulate2(h0, 'b') h1 = -h1[::-1] elif fname == 'test': h0 = np.array([[0, 1, 0], [1, 4, 1], [0, 1, 0]]) h1 = np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]]) elif fname == 'testDVM': h0 = np.array([[1, 1], [1, 1]]) / np.sqrt(2) h1 = np.array([[-1, 1], [1, -1]]) / np.sqrt(2) elif fname == 'qmf': # by Lu, antoniou and xu # ideal response window m = 2 n = 2 w1d = np.kaiser(4*m+1, 2.6) w = np.zeros((n+m+1,n+m+1)) for n1 in np.arange(-m,m+1): for n2 in np.arange(-n,n+1): w[n1+m,n2+n] = w1d[2*m+n1+n2]*w1d[2*m+n1-n2] h = np.zeros((n+m+1,n+m+1)) for n1 in np.arange(-m,m+1): for n2 in np.arange(-n,n+1): h[n1+m, n2+n] = 0.5*np.sinc((n1+n2)/2) * 0.5*np.sinc((n1-n2)/2) c = np.sum(h) h = np.sqrt(2) * h h0 = h * w h1 = modulate2(h0, 'b') elif fname == 'qmf2': # by Lu, Antoniou and Xu # ideal response window h = np.array([ [-0.001104, 0.002494, -0.001744, 0.004895, -0.000048, -0.000311], [0.008918, -0.002844, -0.025197, -0.017135, 0.003905, -0.000081], [-0.007587, -0.065904, 00.100431, -0.055878, 0.007023, 0.001504], [0.001725, 0.184162, 0.632115, 0.099414, -0.027006, -0.001110], [-0.017935, -0.000491, 0.191397, -0.001787, -0.010587, 0.002060], [0.001353, 0.005635, -0.001231, -0.009052, -0.002668, 0.000596]]) h0 = h/np.sum(h) h1 = modulate2(h0, 'b') elif fname == 'dmaxflat4': M1 = 1/np.sqrt(2) M2 = np.copy(M1) k1 = 1-np.sqrt(2) k3 = np.copy(k1) k2 = np.copy(M1) h = np.array([0.25*k2*k3, 0.5*k2, 1+0.5*k2*k3])*M1 h = np.append(h, h[-2::-1]) g = np.array([-0.125*k1*k2*k3, 0.25*k1*k2, -0.5*k1-0.5*k3-0.375*k1*k2*k3, 1+0.5*k1*k2])*M2 g = np.append(g, h[-2::-1]) B = dmaxflat(4,0) h0 = mctrans(h,B) g0 = mctrans(g,B) h0 = np.sqrt(2) * h0 / np.sum(h0) g0 = np.sqrt(2) * g0 / np.sum(g0) h1 = modulate2(g0, 'b') if type == 'r': h1 = modulate2(h0, 'b') h0 = g0 elif fname == 'dmaxflat5': M1 = 1/np.sqrt(2) M2 = M1 k1 = 1-np.sqrt(2) k3 = k1 k2 = M1 h = np.array([0.25*k2*k3, 0.5*k2, 1+0.5*k2*k3])*M1 h = np.append(h, h[-2::-1]) g = np.array([-0.125*k1*k2*k3, 0.25*k1*k2, -0.5*k1-0.5*k3-0.375*k1*k2*k3, 1+0.5*k1*k2])*M2 g = np.append(g, h[-2::-1]) B = dmaxflat(5,0) h0 = mctrans(h,B) g0 = mctrans(g,B) h0 = np.sqrt(2) * h0 / np.sum(h0) g0 = np.sqrt(2) * g0 / np.sum(g0) h1 = modulate2(g0, 'b') if type == 'r': h1 = modulate2(h0, 'b') h0 = g0 elif fname == 'dmaxflat6': M1 = 1/np.sqrt(2) M2 = M1 k1 = 1-np.sqrt(2) k3 = k1 k2 = M1 h = np.array([0.25*k2*k3, 0.5*k2, 1+0.5*k2*k3])*M1 h = np.append(h, h[-2::-1]) g = np.array([-0.125*k1*k2*k3, 0.25*k1*k2, -0.5*k1-0.5*k3-0.375*k1*k2*k3, 1+0.5*k1*k2])*M2 g = np.append(g, h[-2::-1]) B = dmaxflat(6,0) h0 = mctrans(h,B) g0 = mctrans(g,B) h0 = np.sqrt(2) * h0 / np.sum(h0) g0 = np.sqrt(2) * g0 / np.sum(g0) h1 = modulate2(g0, 'b') if type == 'r': h1 = modulate2(h0, 'b') h0 = g0 elif fname == 'dmaxflat7': M1 = 1/np.sqrt(2) M2 = M1 k1 = 1-np.sqrt(2) k3 = k1 k2 = M1 h = np.array([0.25*k2*k3, 0.5*k2, 1+0.5*k2*k3])*M1 h = np.append(h, h[-2::-1]) g = np.array([-0.125*k1*k2*k3, 0.25*k1*k2, -0.5*k1-0.5*k3-0.375*k1*k2*k3, 1+0.5*k1*k2])*M2 g = np.append(g, h[-2::-1]) B = dmaxflat(7,0) h0 = mctrans(h,B) g0 = mctrans(g,B) h0 = np.sqrt(2) * h0 / np.sum(h0) g0 = np.sqrt(2) * g0 / np.sum(g0) h1 = modulate2(g0, 'b') if type == 'r': h1 = modulate2(h0, 'b') h0 = g0 # The original file supports a case "otherwise" for unrecognized filters # and computes simple 1D wavelet filters for them using wfilters.m # I think we don't need this and skip this for the time being. # IN ORIGINAL MATLAB VERSION: # otherwise # % Assume the "degenerated" case: 1D wavelet filters # [h0,h1] = wfilters(fname, type); return h0, h1
def dmaxflat(
N, d)
def dmaxflat(N,d): if (N > 7) or (N < 1): print('Error: N must be in {1,2,...,7}') return 0 if N == 1: h = np.array([[0, 1, 0],[1, 0, 1],[0, 1, 0]])/4 h[1,1] = d elif N == 2: h = np.array([[0, -1, 0],[-1, 0, 10], [0, 10, 0]]) h = np.append(h, np.fliplr(h[:,0:-1]), 1) h = np.append(h, np.flipud(h[0:-1,:]), 0)/32 h[2,2] = d elif N == 3: h = np.array([[0, 3, 0, 2],[3, 0, -27, 0],[0, -27, 0, 174], [2, 0, 174, 0]]) h = np.append(h, np.fliplr(h[:, 0:-1]), 1) h = np.append(h, np.flipud(h[0:-1,:]),0) h[3,3] = d elif N == 4: h = np.array([[0, -5, 0, -3, 0], [-5, 0, 52, 0, 34], [0, 52, 0, -276, 0], [-3, 0, -276, 0, 1454], [0, 34, 0, 1454, 0]])/np.power(2,12) h = np.append(h, np.fliplr(h[:,0:-1]),1) h = np.append(h, np.flipud(h[0:-1,:]),0) h[4,4] = d elif N == 5: h = np.array([[0, 35, 0, 20, 0, 18], [35, 0, -425, 0, -250, 0], [0, -425, 0, 2500, 0, 1610], [20, 0, 2500, 0, -10200, 0], [0, -250, 0, -10200, 0, 47780], [18, 0, 1610, 0, 47780, 0]])/np.power(2,17) h = np.append(h, np.fliplr(h[:,0:-1]),1) h = np.append(h, np.flipud(h[0:-1,:]),0) h[5,5] = d elif N == 6: h = np.array([[0, -63, 0, -35, 0, -30, 0], [-63, 0, 882, 0, 495, 0, 444], [0, 882, 0, -5910, 0, -3420, 0], [-35, 0, -5910, 0, 25875, 0, 16460], [0, 495, 0, 25875, 0, -89730, 0], [-30, 0, -3420, 0, -89730, 0, 389112], [0, 44, 0, 16460, 0, 389112, 0]])/np.power(2,20) h = np.append(h, np.fliplr(h[:,0:-1]),1) h = np.append(h, np.flipud(h[0:-1,:]),0) h[6,6] = d elif N == 7: h = np.array([[0, 231, 0, 126, 0, 105, 0, 100], [231, 0, -3675, 0, -2009, 0, -1715, 0], [0, -3675, 0, 27930, 0, 15435, 0, 13804], [126, 0, 27930, 0, -136514, 0, -77910, 0], [0, -2009, 0, -136514, 0, 495145, 0, 311780], [105, 0, 15435, 0, 495145, 0, -1535709, 0], [0, -1715, 0, -77910, 0, -1534709, 0, 6305740], [100, 0, 13804, 0, 311780, 0, 6305740, 0]])/np.power(2,24) h = np.append(h, np.fliplr(h[:,0:-1]),1) h = np.append(h, np.flipud(h[0:-1,:]),0) h[7,7] = d return h
def mctrans(
b, t)
def mctrans(b,t): # Convert the 1-D filter b to SUM_n a(n) cos(wn) form # if mod(n,2) != 0 -> error n = (b.size-1)//2 b = np.fft.fftshift(b[::-1]) #inverse fftshift b = b[::-1] a = np.zeros(n+1) a[0] = b[0] a[1:n+1] = 2*b[1:n+1] inset = np.floor((np.asarray(t.shape)-1)/2) inset = inset.astype(int) # Use Chebyshev polynomials to compute h P0 = 1 P1 = t; h = a[1]*P1; rows = int(inset[0]+1) cols = int(inset[1]+1) h[rows-1,cols-1] = h[rows-1,cols-1]+a[0]*P0; for i in range(3,n+2): P2 = 2*signal.convolve2d(t, P1) rows = (rows + inset[0]).astype(int) cols = (cols + inset[1]).astype(int) if i == 3: P2[rows-1,cols-1] = P2[rows-1,cols-1] - P0 else: P2[rows[0]-1:rows[-1],cols[0]-1:cols[-1]] = P2[rows[0]-1:rows[-1], cols[0]-1:cols[-1]] - P0 rows = inset[0] + np.arange(np.asarray(P1.shape)[0])+1 rows = rows.astype(int) cols = inset[1] + np.arange(np.asarray(P1.shape)[1])+1 cols = cols.astype(int) hh = h; h = a[i-1]*P2 h[rows[0]-1:rows[-1], cols[0]-1:cols[-1]] = h[rows[0]-1:rows[-1], cols[0]-1:cols[-1]] + hh P0 = P1; P1 = P2; h = np.rot90(h,2) return h
def modulate2(
x, type, center=array([0, 0]))
def modulate2(x, type, center=np.array([0, 0])): size = np.asarray(x.shape) if x.ndim == 1: if np.array_equal(center, [0, 0]): center = 0 origin = np.floor(size/2)+1+center n1 = np.arange(size[0])-origin[0]+1 if x.ndim == 2: n2 = np.arange(size[1])-origin[1]+1 else: n2 = n1 if type == 'r': m1 = np.power(-1,n1) if x.ndim == 1: y = x*m1 else: y = x * np.transpose(np.tile(m1, (size[1], 1))) elif type == 'c': m2 = np.power(-1,n2) if x.ndim == 1: y = x*m2 else: y = x * np.tile(m2, np.array([size[0], 1])) elif type == 'b': m1 = np.power(-1,n1) m2 = np.power(-1,n2) m = np.outer(m1, m2) if x.ndim == 1: y = x * m1 else: y = x * m return y