Wavelets and polyphase transfer / exact Python program
Four Haar levels, circular three-tap FIR filtering, exact polyphase and all-four-level Haar operator transfers, exact lifting, weighted coefficient budgets, finite filter stability, order-of-operations comparison, boundary repair and counterexamples. Python 3.10 or later, standard library only.
Read the complete worked example Download Python · Source and results bundle · Source identities
Reproducible records: phase-transfer.csv · haar-coefficients.csv · weighted-budget.csv · full-wavelet-transfer.csv · multiscale-energy.csv · operation-order.csv · results.json · signal-comparison.svg · full-wavelet-transfer.svg
Run python -B wavelet_polyphase.py. Add --output wavelet-results to generate all exact CSV, SVG, and JSON records.
1 #!/usr/bin/env python3
2 """Exact finite Haar and polyphase FIR transfer, with reproducible source data.
3
4 All arithmetic uses Fraction. The FIR filter is circular on a declared finite
5 window; no infinite-signal or zero-extension claim is made.
6 """
7 from __future__ import annotations
8
9 import argparse
10 import csv
11 from fractions import Fraction as F
12 from itertools import product
13 import json
14 from pathlib import Path
15 from random import Random
16
17 DEFAULT = (2, 2, 2, 2, 6, 6, 6, 6, -4, -4, 12, -4, 4, 4, 4, 4)
18 TAPS = (1, 2, -1)
19
20
21 def exact(values):
22 return tuple(F(v) for v in values)
23
24
25 def split_phases(values, modulus=2):
26 values = exact(values)
27 if modulus < 1 or not values or len(values) % modulus:
28 raise ValueError("The nonempty signal length must be divisible by the phase count")
29 return tuple(values[r::modulus] for r in range(modulus))
30
31
32 def join_phases(phases):
33 phases = tuple(exact(row) for row in phases)
34 if not phases or not phases[0] or any(len(row) != len(phases[0]) for row in phases):
35 raise ValueError("Phase rows must be nonempty and equally long")
36 return tuple(phases[r][k] for k in range(len(phases[0])) for r in range(len(phases)))
37
38
39 def shift(values):
40 """The periodic right shift (T x)[i] = x[(i-1) mod len(x)]."""
41 values = exact(values)
42 if not values:
43 raise ValueError("Signal must be nonempty")
44 return values[-1:] + values[:-1]
45
46
47 def backward_shift(values):
48 """The inverse cyclic shift (T^-1 x)[i] = x[(i+1) mod len(x)]."""
49 values = exact(values)
50 if not values:
51 raise ValueError("Signal must be nonempty")
52 return values[1:] + values[:1]
53
54
55 def skew_shift(values):
56 """(T^-1 - T)x, the finite energy certificate for the FIR."""
57 return tuple(left - right for left, right
58 in zip(backward_shift(values), shift(values)))
59
60
61 def zero_extended_fir(values, taps=TAPS):
62 """Finite causal convolution with values outside the window set to zero."""
63 values = exact(values)
64 if not values:
65 raise ValueError("Signal must be nonempty")
66 return tuple(sum((F(tap) * values[i-j]
67 for j, tap in enumerate(taps) if 0 <= i-j < len(values)),
68 F(0)) for i in range(len(values)))
69
70
71 def fir(values, taps=TAPS):
72 """Circular FIR: H = sum_j taps[j] T^j."""
73 values = exact(values)
74 if not values:
75 raise ValueError("Signal must be nonempty")
76 out = (F(0),) * len(values)
77 shifted = values
78 for coefficient in taps:
79 out = tuple(a + F(coefficient) * b for a, b in zip(out, shifted))
80 shifted = shift(shifted)
81 return out
82
83
84 def haar_pair(even, odd):
85 even, odd = exact(even), exact(odd)
86 if len(even) != len(odd):
87 raise ValueError("Haar phases need equal length")
88 return (tuple((u + v) / 2 for u, v in zip(even, odd)),
89 tuple((u - v) / 2 for u, v in zip(even, odd)))
90
91
92 def lifting_decompose(even, odd):
93 """Paper V predict/update: difference=odd-even, average=even+difference/2."""
94 even, odd = exact(even), exact(odd)
95 if len(even) != len(odd):
96 raise ValueError("Lifting phases need equal length")
97 differences = tuple(o-e for e, o in zip(even, odd))
98 averages = tuple(e + difference/2
99 for e, difference in zip(even, differences))
100 return averages, differences
101
102
103 def lifting_synthesize(averages, differences):
104 """Invert the predict/update map using exactly the stored side channel."""
105 averages, differences = exact(averages), exact(differences)
106 if len(averages) != len(differences):
107 raise ValueError("Lifting phases need equal length")
108 even = tuple(a - difference/2
109 for a, difference in zip(averages, differences))
110 odd = tuple(e + difference for e, difference in zip(even, differences))
111 return join_phases((even, odd))
112
113
114 def inverse_haar_pair(averages, details):
115 averages, details = exact(averages), exact(details)
116 if len(averages) != len(details):
117 raise ValueError("Haar phases need equal length")
118 return join_phases((tuple(a + d for a, d in zip(averages, details)),
119 tuple(a - d for a, d in zip(averages, details))))
120
121
122 def analyze(values):
123 """Return final mean and all (averages, details), fine to coarse."""
124 values = exact(values)
125 if len(values) < 2 or len(values) & (len(values) - 1):
126 raise ValueError("Haar analysis needs a power-of-two length >= 2")
127 current, levels = values, []
128 while len(current) > 1:
129 averages, details = haar_pair(*split_phases(current))
130 levels.append((averages, details))
131 current = averages
132 return current[0], tuple(levels)
133
134
135 def synthesize(mean, levels, omit=frozenset()):
136 """Reconstruct. An omitted detail (level, index) is replaced by zero."""
137 current = (F(mean),)
138 for j in range(len(levels), 0, -1):
139 details = levels[j - 1][1]
140 if len(current) != len(details):
141 raise ValueError("Inconsistent Haar level lengths")
142 saved = tuple(F(0) if (j, k) in omit else d for k, d in enumerate(details))
143 current = inverse_haar_pair(current, saved)
144 return current
145
146
147 def energy(values):
148 return sum((F(v) ** 2 for v in values), F(0))
149
150
151 def coefficient_energy(mean, levels):
152 length = 2 ** len(levels)
153 return length * F(mean) ** 2 + sum(
154 (2 ** j * energy(pair[1]) for j, pair in enumerate(levels, 1)), F(0))
155
156
157 def best_detail_budget(levels):
158 """Optimal fixed Haar-coordinate budget with the mean always retained.
159
160 The detail synthesis vectors are orthogonal with squared norm 2**level.
161 Sort exact weighted energies; ties use (level,index) to be deterministic.
162 Returns (retained_count, retained_energy, residual_energy) for each budget.
163 This is not an adaptive best-packet-tree algorithm.
164 """
165 entries = sorted(
166 ((2**j * F(d)**2, j, k)
167 for j, (_, details) in enumerate(levels, 1)
168 for k, d in enumerate(details) if d != 0),
169 key=lambda item: (-item[0], item[1], item[2]))
170 residual = sum((item[0] for item in entries), F(0))
171 records = [(0, F(0), residual)]
172 retained = F(0)
173 for budget, (contribution, _, _) in enumerate(entries, 1):
174 retained += contribution
175 residual -= contribution
176 records.append((budget, retained, residual))
177 return tuple(records)
178
179
180 def wavelet_analysis(values):
181 """All lowpass-branch Haar coordinates: mean, then coarse-to-fine details."""
182 mean, levels = analyze(values)
183 return (mean,) + tuple(d for _, details in reversed(levels) for d in details)
184
185
186 def wavelet_synthesis(coordinates):
187 """Inverse transform using only the declared coarse-to-fine coordinates."""
188 coordinates = exact(coordinates)
189 n = len(coordinates)
190 if n < 2 or n & (n-1):
191 raise ValueError("Wavelet coordinates need a power-of-two length >= 2")
192 current = (coordinates[0],)
193 offset = 1
194 while len(current) < n:
195 count = len(current)
196 current = inverse_haar_pair(current, coordinates[offset:offset+count])
197 offset += count
198 if offset != n:
199 raise AssertionError("Wavelet coordinate partition did not close")
200 return current
201
202
203 def wavelet_labels(length):
204 if length < 2 or length & (length - 1):
205 raise ValueError("Wavelet coordinate labels require a dyadic length")
206 depth = length.bit_length()-1
207 return ("m",) + tuple(f"L{j}.{k}" for j in range(depth,0,-1)
208 for k in range(length//(2**j)))
209
210
211 def wavelet_energy_levels(coordinates):
212 """Five weighted bands for length 16; generalized to any dyadic length."""
213 coordinates = exact(coordinates)
214 labels = wavelet_labels(len(coordinates))
215 depth = len(coordinates).bit_length()-1
216 result = [("mean", len(coordinates)*coordinates[0]**2,
217 int(coordinates[0] != 0))]
218 offset = 1
219 for j in range(depth,0,-1):
220 width = len(coordinates)//(2**j)
221 values = coordinates[offset:offset+width]
222 result.append((f"level_{j}", (2**j)*energy(values),
223 sum(int(value != 0) for value in values)))
224 offset += width
225 if len(labels) != offset:
226 raise AssertionError("The scale labels do not match the input")
227 return tuple(result)
228
229
230 def full_wavelet_transfer(length=16):
231 """Exact W H W^-1 matrix, columns obtained from wavelet coordinate impulses."""
232 wavelet_labels(length)
233 columns = []
234 for j in range(length):
235 e = tuple(F(k == j) for k in range(length))
236 columns.append(wavelet_analysis(fir(wavelet_synthesis(e))))
237 return tuple(tuple(columns[j][i] for j in range(length))
238 for i in range(length))
239
240
241 def apply_matrix(matrix, values):
242 values = exact(values)
243 if len(matrix) != len(values) or any(len(row) != len(values) for row in matrix):
244 raise ValueError("Transfer matrix and coordinate length differ")
245 return tuple(sum((F(a)*b for a,b in zip(row,values)),F(0)) for row in matrix)
246
247
248 def coarsen_first_level(values):
249 """Project onto first-level Haar averages, duplicating every pair mean."""
250 values = exact(values)
251 if len(values)<2 or len(values)%2:
252 raise ValueError("First-level coarsening needs a positive even length")
253 return inverse_haar_pair(haar_pair(*split_phases(values))[0],
254 (F(0),)*(len(values)//2))
255
256
257 def operation_commutator(values):
258 """(P H - H P)x, where P is the averaging projection."""
259 values = exact(values)
260 ph = coarsen_first_level(fir(values))
261 hp = fir(coarsen_first_level(values))
262 return tuple(a-b for a,b in zip(ph,hp))
263
264
265 def multiscale_transfer_svg(matrix):
266 """Python-generated signed rational heatmap with matching coordinate labels."""
267 n = len(matrix)
268 labels = wavelet_labels(n)
269 if len(matrix)!=n or any(len(row)!=n for row in matrix):
270 raise ValueError("Heatmap requires a square wavelet matrix")
271 cell = 34
272 left, top = 138, 112
273 width, height = 822, 774
274 nnz=sum(int(coefficient != 0) for row in matrix for coefficient in row)
275 parts = [
276 '<svg xmlns="http://www.w3.org/2000/svg" width="822" height="774" viewBox="0 0 822 774" role="img" aria-labelledby="title description">',
277 '<title id="title">Exact four-level Haar transfer matrix</title>',
278 '<desc id="description">A 16 by 16 matrix for W H W inverse. Rows are output wavelet coordinates, columns are input wavelet coordinates, ordered global mean then levels four through one. Teal is positive, ochre negative, grey zero. Each nonzero cell is labeled with its exact rational value.</desc>',
279 '<rect width="822" height="774" fill="#fff"/>',
280 '<text x="36" y="34" font-family="sans-serif" font-size="21" font-weight="600" fill="#173f3a">FIR transfer in four-level Haar coordinates</text>',
281 f'<text x="36" y="60" font-family="sans-serif" font-size="14" fill="#425852">{n} × {n} exact coefficients · {nnz} nonzero · W H W⁻¹</text>',
282 '<text x="140" y="85" font-family="sans-serif" font-size="12" fill="#536b66">Input wavelet coordinates (columns)</text>',
283 '<text x="12" y="107" font-family="sans-serif" font-size="12" fill="#536b66">Output rows</text>'
284 ]
285 for j, label in enumerate(labels):
286 cx = left + j*cell + cell//2
287 parts.append(f'<text x="{cx}" y="104" text-anchor="middle" font-family="monospace" font-size="10" fill="#445f59">{label}</text>')
288 for i, row in enumerate(matrix):
289 y = top+i*cell
290 parts.append(f'<text x="{left-10}" y="{y+22}" text-anchor="end" font-family="monospace" font-size="11" fill="#334e49">{labels[i]}</text>')
291 for j, value in enumerate(row):
292 x = left+j*cell
293 fill="#edf1f0" if value==0 else ("#28786d" if value>0 else "#965023")
294 parts.append(f'<rect x="{x}" y="{y}" width="{cell-1}" height="{cell-1}" rx="2" fill="{fill}"/>')
295 if value:
296 foreground="#ffffff"
297 parts.append(f'<text x="{x+cell//2}" y="{y+21}" text-anchor="middle" font-family="monospace" font-size="10" font-weight="600" fill="{foreground}">{value}</text>')
298 parts.extend((
299 '<rect x="140" y="686" width="16" height="16" fill="#28786d"/>',
300 '<text x="165" y="699" font-family="sans-serif" font-size="13" fill="#334e49">positive</text>',
301 '<rect x="258" y="686" width="16" height="16" fill="#965023"/>',
302 '<text x="283" y="699" font-family="sans-serif" font-size="13" fill="#334e49">negative</text>',
303 '<rect x="382" y="686" width="16" height="16" fill="#edf1f0"/>',
304 '<text x="407" y="699" font-family="sans-serif" font-size="13" fill="#334e49">zero</text>',
305 '<text x="140" y="740" font-family="sans-serif" font-size="13" fill="#536b66">Averaging-normalized Haar basis · no approximation · finite circular boundary</text>',
306 '</svg>'
307 ))
308 return "\n".join(parts) + "\n"
309
310
311 def transfer_polyphase(even, odd):
312 """(E',O') = ((I-Z)E+2ZO, 2E+(I-Z)O), Z a circular half-rate shift."""
313 even, odd = exact(even), exact(odd)
314 if len(even) != len(odd) or not even:
315 raise ValueError("Equal nonempty phase rows required")
316 ze, zo = shift(even), shift(odd)
317 return (tuple(e - ep + 2 * op for e, ep, op in zip(even, ze, zo)),
318 tuple(2 * e + o - op for e, o, op in zip(even, odd, zo)))
319
320
321 def transfer_haar(averages, details):
322 """(A',D') = (2A+(I-Z)D, (Z-I)A-2ZD)."""
323 averages, details = exact(averages), exact(details)
324 if len(averages) != len(details) or not averages:
325 raise ValueError("Equal nonempty coefficient rows required")
326 za, zd = shift(averages), shift(details)
327 return (tuple(2*a + d - prev for a, d, prev in zip(averages, details, zd)),
328 tuple(prev - a - 2*old for a, prev, old in zip(averages, za, zd)))
329
330
331 def diagonal_only(averages, details):
332 """Deliberately invalid shortcut: erase the two off-diagonal transfers."""
333 a, d = exact(averages), exact(details)
334 return inverse_haar_pair(tuple(2*x for x in a), tuple(-2*x for x in shift(d)))
335
336
337 def default_case():
338 x = exact(DEFAULT)
339 mean, levels = analyze(x)
340 omit = frozenset((1, k) for k in range(len(levels[0][1])))
341 coarse = synthesize(mean, levels, omit)
342 residual = tuple(u - v for u, v in zip(x, coarse))
343 filtered = fir(x)
344 even, odd = split_phases(x)
345 a, d = haar_pair(even, odd)
346 wrong = diagonal_only(a, d)
347 boundary = zero_extended_fir(x)
348 return dict(x=x, mean=mean, levels=levels, coarse=coarse, residual=residual,
349 filtered=filtered, wrong=wrong, boundary=boundary,
350 boundary_residual=tuple(u - v for u, v in zip(filtered, boundary)),
351 filter_error=tuple(u - v for u, v in zip(filtered, wrong)))
352
353
354 def check():
355 counts = {}
356 def require(condition, group, message):
357 if not condition:
358 raise AssertionError(group + ": " + message)
359 counts[group] = counts.get(group, 0) + 1
360
361 case = default_case()
362 x, mean, levels = case["x"], case["mean"], case["levels"]
363 require(mean == 3, "reference", "global mean")
364 for j, expected in enumerate((
365 (0, 0, 0, 0, 0, 8, 0, 0),
366 (0, 0, -4, 0),
367 (-2, -2),
368 (1,)), 1):
369 require(levels[j-1][1] == exact(expected), "reference", "detail level "+str(j))
370 require(synthesize(mean, levels) == x, "reference", "full inverse")
371 require(energy(x) == coefficient_energy(mean, levels) == 416,
372 "reference", "weighted energy")
373 require(case["coarse"] == exact((2,2,2,2,6,6,6,6,-4,-4,4,4,4,4,4,4)),
374 "reference", "coarse sample sequence")
375 require(case["residual"] == exact((0,0,0,0,0,0,0,0,0,0,8,-8,0,0,0,0)),
376 "reference", "missing spike detail")
377 require(energy(case["coarse"]) == 288 and energy(case["residual"]) == 128
378 and energy(x) == energy(case["coarse"]) + energy(case["residual"]),
379 "reference", "exact Pythagorean split")
380 require(case["filtered"] == exact((6,2,4,4,8,16,12,12,2,-18,8,24,-16,16,8,8)),
381 "reference", "circular FIR samples")
382 require(energy(case["filtered"]) == 2288, "reference", "filtered energy")
383 require(energy(case["filter_error"]) == 624, "reference", "off-diagonal omission")
384 require(fir(x)[0] == x[0] + 2*x[-1] - x[-2] == 6,
385 "reference", "periodic boundary")
386 require(split_phases(x, 4) == (
387 exact((2,6,-4,4)),exact((2,6,-4,4)),
388 exact((2,6,12,4)),exact((2,6,-4,4))),
389 "reference", "four-address refinement")
390 require(join_phases(split_phases(x, 2)) == join_phases(split_phases(x, 4)) == x,
391 "reference", "two and four address synthesis")
392 require(haar_pair(*split_phases(case["filtered"])) ==
393 transfer_haar(*haar_pair(*split_phases(x))),
394 "reference", "transferred Haar coefficients")
395 require(join_phases(transfer_polyphase(*split_phases(x))) == case["filtered"],
396 "reference", "polyphase transfer")
397 require(diagonal_only(*haar_pair(*split_phases(x))) != case["filtered"],
398 "obstructions", "cross phase interactions cannot be discarded")
399 require(haar_pair(*split_phases(case["coarse"]))[0] ==
400 haar_pair(*split_phases(x))[0], "obstructions", "lowpass collision")
401 require(case["coarse"] != x, "obstructions", "different signals same lowpass")
402 require(synthesize(*analyze(case["coarse"])) == case["coarse"],
403 "obstructions", "coarsened signal is valid")
404 even, odd = split_phases(x)
405 require(join_phases((even, (F(0),)*len(odd))) != x,
406 "obstructions", "missing odd phase")
407 require(x[10] == 12 and x[11] == -4 and
408 case["filtered"][10] == 8 and case["filtered"][11] == 24,
409 "reference", "local boundary-crossing calculation")
410
411 # The independent Paper V lifting factorization uses d_L = odd-even.
412 a, d = haar_pair(*split_phases(x))
413 lifted_a, lifted_difference = lifting_decompose(*split_phases(x))
414 require(lifted_a == a and lifted_difference == tuple(-2*v for v in d),
415 "reference", "lifting equals Haar after declared rescaling")
416 require(lifting_synthesize(lifted_a, lifted_difference) == x,
417 "reference", "exact predict/update synthesis")
418 skew = skew_shift(x)
419 require(energy(skew) == 624 and
420 energy(case["filtered"]) == 4*energy(x) + energy(skew),
421 "reference", "sharp FIR energy identity")
422 require(4*energy(x) <= energy(case["filtered"]) <= 8*energy(x),
423 "reference", "finite stability interval")
424 require(energy(case["filter_error"]) == energy(skew),
425 "reference", "off-diagonal error equals FIR energy surplus")
426 require(case["boundary_residual"] == exact((4,-4)+(0,)*14)
427 and energy(case["boundary_residual"]) == 32,
428 "reference", "two-site boundary correction")
429 require(tuple(record[2] for record in best_detail_budget(levels))
430 == tuple(map(F, (272,144,80,48,16,0))),
431 "reference", "optimal fixed-basis detail budget")
432 contrast = exact((2,)*8 + (-2,-2,1,-5,-2,-2,-2,-2))
433 contrast_mean, contrast_levels = analyze(contrast)
434 require(contrast_mean == 0 and energy(contrast) == 82
435 and best_detail_budget(contrast_levels)[1] == (1,F(64),F(18)),
436 "reference", "weighted selection beats largest raw magnitude")
437 flat, quarter = exact((1,)*16), exact((1,0,-1,0)*4)
438 require(energy(fir(flat)) == 4*energy(flat) and
439 energy(fir(quarter)) == 8*energy(quarter),
440 "reference", "both finite operator norm bounds are sharp")
441
442 # Full four-level operator transfer, not merely one Haar pair.
443 wx = wavelet_analysis(x)
444 wy = wavelet_analysis(case["filtered"])
445 matrix16 = full_wavelet_transfer(16)
446 project = coarsen_first_level
447 ph = project(case["filtered"])
448 hp = fir(project(x))
449 commutator = operation_commutator(x)
450 require(wx == exact((3,1,-2,-2,0,0,-4,0,0,0,0,0,0,8,0,0)),
451 "multiscale", "full input coefficient ordering")
452 require(wy == exact((6,2,-4,0,0,0,-12,-4,2,0,-4,0,10,-8,-16,0)),
453 "multiscale", "full filtered coefficient ordering")
454 require(wavelet_synthesis(wx) == x,
455 "multiscale", "wavelet analysis-synthesis inverse")
456 require(wavelet_synthesis(wy) == case["filtered"],
457 "multiscale", "filtered wavelet analysis-synthesis")
458 require(tuple(value for _,value,_ in wavelet_energy_levels(wx))
459 == exact((144,16,64,64,128)),
460 "multiscale", "five original weighted energy bands")
461 require(tuple(value for _,value,_ in wavelet_energy_levels(wy))
462 == exact((576,64,128,640,880)),
463 "multiscale", "five filtered weighted energy bands")
464 require(sum(int(t!=0) for t in wx) == 6,
465 "multiscale", "original wavelet support size")
466 require(sum(int(t!=0) for t in wy) == 10,
467 "multiscale", "filtered wavelet support size")
468 require(len(matrix16) == 16 and all(len(row)==16 for row in matrix16),
469 "multiscale", "full square transfer matrix shape")
470 require(sum(t!=0 for row in matrix16 for t in row) == 56,
471 "multiscale", "exact transfer matrix sparsity")
472 require(apply_matrix(matrix16, wx) == wy,
473 "multiscale", "full four-level operator transfer on the signal")
474 require(matrix16[0] == exact((2,)+(0,)*15),
475 "multiscale", "the global mean scales exactly by two")
476 require(project(x) == case["coarse"],
477 "multiscale", "coarsening equals dropping finest details")
478 require(project(project(x)) == project(x),
479 "multiscale", "coarsening is idempotent")
480 require(energy(commutator) == 624,
481 "multiscale", "noncommutation witness energy")
482 require(energy(tuple(case["filtered"][i]-hp[i] for i in range(16))) == 768
483 and energy(tuple(case["filtered"][i]-ph[i] for i in range(16))) == 880,
484 "multiscale", "distinct orders have distinct projection errors")
485 require(energy(commutator) == energy(case["filter_error"]),
486 "multiscale", "commutator norm equals omitted-transfer energy")
487 require(commutator == exact((-2,2,0,0,4,-4,0,0,-10,10,16,0,-8,-8,0,0)),
488 "multiscale", "exact commutator coordinate witness")
489 require(commutator != (F(0),)*16,
490 "multiscale", "filter and coarse projection do not commute")
491
492 # Linearity plus equality on every coordinate basis vector gives an
493 # exact 16-dimensional operator identity, not merely a sample agreement.
494 for k in range(16):
495 basis = tuple(F(i == k) for i in range(16))
496 e, o = split_phases(basis)
497 a, d = haar_pair(e, o)
498 direct = fir(basis)
499 require(join_phases(transfer_polyphase(e, o)) == direct,
500 "basis", "polyphase FIR basis "+str(k))
501 require(haar_pair(*split_phases(direct)) == transfer_haar(a, d),
502 "basis", "Haar conjugacy basis "+str(k))
503 require(synthesize(*analyze(basis)) == basis, "basis", "wavelet basis "+str(k))
504 lifted_a, lifted_d = lifting_decompose(e, o)
505 require(lifted_a == a and lifted_d == tuple(-2*v for v in d)
506 and lifting_synthesize(lifted_a, lifted_d) == basis,
507 "basis", "lifting operator basis "+str(k))
508 require(energy(direct) == 4*energy(basis)+energy(skew_shift(basis)),
509 "basis", "FIR adjoint energy basis "+str(k))
510 require(energy(tuple(u-v for u,v in zip(direct, diagonal_only(a,d))))
511 == energy(skew_shift(basis)),
512 "basis", "omitted cross phase error basis "+str(k))
513 residue = tuple(u-v for u,v in zip(direct, zero_extended_fir(basis)))
514 require(residue == (2*basis[-1]-basis[-2], -basis[-1])
515 + (F(0),)*(len(basis)-2),
516 "basis", "boundary correction basis "+str(k))
517 waves = wavelet_analysis(basis)
518 require(wavelet_synthesis(waves) == basis,
519 "multiscale", "four-level inverse on coordinate impulse")
520 require(apply_matrix(matrix16, waves) == wavelet_analysis(direct),
521 "multiscale", "whole wavelet matrix on coordinate impulse")
522 require(energy(operation_commutator(basis))
523 == energy(skew_shift(basis)),
524 "multiscale", "whole commutator on sample impulse")
525 require(energy(basis) == sum(v for _,v,_ in wavelet_energy_levels(waves)),
526 "multiscale", "wavelet energy on sample impulse")
527
528 # Exhaustive small-integer 4-vectors and deterministic longer controls.
529 matrix4=full_wavelet_transfer(4)
530 for numbers in product((-2, -1, 0, 1, 2), repeat=4):
531 s = exact(numbers)
532 m, stages = analyze(s)
533 require(synthesize(m, stages) == s, "exhaustive", "four-sample inverse")
534 require(coefficient_energy(m, stages) == energy(s),
535 "exhaustive", "four-sample energy")
536 require(join_phases(transfer_polyphase(*split_phases(s))) == fir(s),
537 "exhaustive", "four-sample polyphase FIR")
538 even, odd = split_phases(s)
539 a, d = haar_pair(even, odd)
540 lifted_a, lifted_d = lifting_decompose(even, odd)
541 require(lifted_a == a and lifted_d == tuple(-2*v for v in d)
542 and lifting_synthesize(lifted_a, lifted_d) == s,
543 "exhaustive", "four-sample lifting inverse")
544 direct = fir(s)
545 require(energy(direct) == 4*energy(s) + energy(skew_shift(s)),
546 "exhaustive", "four-sample FIR energy")
547 require(energy(tuple(u-v for u,v in zip(direct, diagonal_only(a,d))))
548 == energy(skew_shift(s)),
549 "exhaustive", "four-sample omitted cross blocks")
550 boundary_difference = tuple(u-v for u,v in zip(direct, zero_extended_fir(s)))
551 require(boundary_difference == (2*s[-1]-s[-2], -s[-1])
552 + (F(0),)*(len(s)-2),
553 "exhaustive", "four-sample edge correction")
554 require(apply_matrix(matrix4, wavelet_analysis(s)) == wavelet_analysis(direct),
555 "multiscale", "exhaustive 4-point whole wavelet transfer")
556 require(energy(operation_commutator(s)) == energy(skew_shift(s)),
557 "multiscale", "exhaustive 4-point noncommutation certificate")
558 rng = Random(20261030)
559 wavelet_matrices = {8:full_wavelet_transfer(8),16:matrix16,32:full_wavelet_transfer(32)}
560 for length in (8, 16, 32):
561 for _ in range(32):
562 s = exact(rng.randrange(-12, 13) for _ in range(length))
563 m, stages = analyze(s)
564 require(synthesize(m, stages) == s, "random", "multilevel synthesis")
565 require(coefficient_energy(m, stages) == energy(s),
566 "random", "weighted Parseval identity")
567 require(join_phases(split_phases(s, 2)) == s, "random", "phase inverse")
568 require(join_phases(transfer_polyphase(*split_phases(s))) == fir(s),
569 "random", "circular FIR polyphase transfer")
570 require(haar_pair(*split_phases(fir(s))) ==
571 transfer_haar(*haar_pair(*split_phases(s))),
572 "random", "Haar coordinate transfer")
573 even, odd = split_phases(s)
574 a, d = haar_pair(even, odd)
575 lifted_a, lifted_d = lifting_decompose(even, odd)
576 require(lifted_a == a and lifted_d == tuple(-2*v for v in d)
577 and lifting_synthesize(lifted_a, lifted_d) == s,
578 "random", "predict update lifting")
579 direct = fir(s)
580 require(energy(direct) == 4*energy(s) + energy(skew_shift(s)),
581 "random", "sharp finite filter energy")
582 require(energy(tuple(u-v for u,v in zip(direct, diagonal_only(a,d))))
583 == energy(skew_shift(s)),
584 "random", "cross channel energy identity")
585 boundary_difference = tuple(u-v for u,v in zip(direct, zero_extended_fir(s)))
586 require(boundary_difference == (2*s[-1]-s[-2], -s[-1])
587 + (F(0),)*(len(s)-2),
588 "random", "cyclic versus zero edge packet")
589 require(apply_matrix(wavelet_matrices[length], wavelet_analysis(s))
590 == wavelet_analysis(direct),
591 "multiscale", "whole wavelet transfer for seeded signal")
592 require(energy(operation_commutator(s)) == energy(skew_shift(s)),
593 "multiscale", "noncommutation for seeded signal")
594 return counts
595
596
597 def figure_svg(case):
598 """A small, accessible SVG generated from exact sample data."""
599 width, height = 1020, 520
600 colors = ("#205e54", "#287d89", "#b96b35")
601 panels = ((case["x"], "Original x", colors[0]),
602 (case["coarse"], "Coarse reconstruction (finest details removed)", colors[1]),
603 (case["residual"], "Residual x - reconstruction", colors[2]))
604 pieces = [
605 '<svg xmlns="http://www.w3.org/2000/svg" width="1020" height="520" viewBox="0 0 1020 520" role="img" aria-labelledby="title description">',
606 '<title id="title">Haar reconstruction and exact residual</title>',
607 '<desc id="description">Three aligned sixteen-sample stem plots: the input signal, its lowpass reconstruction, and the missing finest detail. The residual has values eight and negative eight only at positions ten and eleven.</desc>',
608 '<rect width="1020" height="520" fill="#ffffff"/>',
609 '<text x="35" y="30" font-family="sans-serif" font-size="20" font-weight="600" fill="#193b36">One signal · full, coarse and residual coordinates</text>',
610 ]
611 def fmt(v):
612 return format(float(v), ".2f")
613 for index, (data, name, color) in enumerate(panels):
614 top = 48 + index * 146
615 baseline = top + 88
616 pieces.extend((
617 f'<text x="36" y="{top+15}" font-family="sans-serif" font-size="15" font-weight="600" fill="#203c39">{name}</text>',
618 f'<line x1="41" y1="{baseline}" x2="972" y2="{baseline}" stroke="#9eb3af" stroke-width="1.2"/>',
619 f'<rect x="30" y="{top+23}" width="950" height="114" rx="6" fill="none" stroke="#e3eeeb"/>',
620 ))
621 for k, value in enumerate(data):
622 px = 60 + 58*k
623 py = baseline - 4.8*float(value)
624 pieces.append(f'<line x1="{px}" y1="{baseline}" x2="{px}" y2="{fmt(py)}" stroke="{color}" stroke-width="3"/>')
625 pieces.append(f'<circle cx="{px}" cy="{fmt(py)}" r="4" fill="{color}"/>')
626 if k in (0, 4, 8, 12, 15):
627 pieces.append(f'<text x="{px}" y="{top+130}" text-anchor="middle" font-family="sans-serif" font-size="12" fill="#526966">{k}</text>')
628 pieces.extend((
629 '<text x="977" y="504" text-anchor="end" font-family="sans-serif" font-size="11" fill="#526966">Sample addresses 0–15 · linear amplitude scale in every panel</text>',
630 '</svg>'))
631 return "\n".join(pieces) + "\n"
632
633
634 def save(folder, counts):
635 folder.mkdir(parents=True, exist_ok=True)
636 case = default_case()
637 x, levels = case["x"], case["levels"]
638 with (folder/"phase-transfer.csv").open("w", encoding="utf-8", newline="") as stream:
639 w = csv.writer(stream, lineterminator="\n")
640 w.writerow(("index", "original", "coarse", "residual", "circular_fir",
641 "incorrect_diagonal_only", "transfer_error",
642 "zero_extended_fir", "cyclic_minus_zero_boundary"))
643 for i in range(len(x)):
644 w.writerow((i, *(str(v[i]) for v in (
645 x, case["coarse"], case["residual"], case["filtered"],
646 case["wrong"], case["filter_error"], case["boundary"],
647 case["boundary_residual"]))))
648 with (folder/"haar-coefficients.csv").open("w", encoding="utf-8", newline="") as stream:
649 w = csv.writer(stream, lineterminator="\n")
650 w.writerow(("kind", "level", "index", "value", "energy_weight", "weighted_energy"))
651 w.writerow(("mean", 4, 0, str(case["mean"]), 16, str(16*case["mean"]**2)))
652 for j in range(len(levels), 0, -1):
653 for k, d in enumerate(levels[j-1][1]):
654 w.writerow(("detail", j, k, str(d), 2**j, str(2**j*d*d)))
655 with (folder/"weighted-budget.csv").open("w", encoding="utf-8", newline="") as stream:
656 w = csv.writer(stream, lineterminator="\n")
657 w.writerow(("retained_details", "retained_detail_energy", "squared_residual_error"))
658 for budget, retained, residual in best_detail_budget(levels):
659 w.writerow((budget, str(retained), str(residual)))
660 input_waves = wavelet_analysis(x)
661 output_waves = wavelet_analysis(case["filtered"])
662 matrix = full_wavelet_transfer(len(x))
663 labels = wavelet_labels(len(x))
664 with (folder/"full-wavelet-transfer.csv").open("w", encoding="utf-8", newline="") as stream:
665 w = csv.writer(stream, lineterminator="\n")
666 w.writerow(("output_index","output_coordinate","input_index","input_coordinate","coefficient"))
667 for i, row in enumerate(matrix):
668 for j, coefficient in enumerate(row):
669 if coefficient:
670 w.writerow((i, labels[i], j, labels[j], str(coefficient)))
671 input_bands = wavelet_energy_levels(input_waves)
672 output_bands = wavelet_energy_levels(output_waves)
673 with (folder/"multiscale-energy.csv").open("w", encoding="utf-8", newline="") as stream:
674 w = csv.writer(stream, lineterminator="\n")
675 w.writerow(("band","original_energy","filtered_energy","original_nonzero","filtered_nonzero"))
676 for (name, en0, count0), (_, en1, count1) in zip(input_bands,output_bands):
677 w.writerow((name, str(en0), str(en1), count0, count1))
678 ph=coarsen_first_level(case["filtered"])
679 hp=fir(coarsen_first_level(x))
680 with (folder/"operation-order.csv").open("w", encoding="utf-8", newline="") as stream:
681 w=csv.writer(stream, lineterminator="\n")
682 w.writerow(("index","original","filtered","coarse_original","filter_coarse",
683 "coarse_filtered","commutator","filtered_minus_filter_coarse",
684 "filtered_minus_coarse_filtered"))
685 for k in range(len(x)):
686 w.writerow((k,str(x[k]),str(case["filtered"][k]),str(case["coarse"][k]),
687 str(hp[k]),str(ph[k]),str(ph[k]-hp[k]),
688 str(case["filtered"][k]-hp[k]),
689 str(case["filtered"][k]-ph[k])))
690 (folder/"full-wavelet-transfer.svg").write_text(
691 multiscale_transfer_svg(matrix), encoding="utf-8", newline="\n")
692 (folder/"signal-comparison.svg").write_text(figure_svg(case), encoding="utf-8", newline="\n")
693 record = {
694 "case": "Haar averaging / cyclic FIR (1, 2, -1)",
695 "signal": list(DEFAULT),
696 "filtered": [int(v) for v in case["filtered"]],
697 "coarse": [int(v) for v in case["coarse"]],
698 "energies": {"original": 416, "coarse": 288, "residual": 128,
699 "filtered": 2288, "missing_transfer_blocks": 624},
700 "stability": {
701 "fourfold_input": 1664,
702 "skew_shift_energy": int(energy(skew_shift(x))),
703 "filtered": int(energy(case["filtered"])),
704 "eightfold_input": 3328,
705 "omitted_transfer_error": int(energy(case["filter_error"])),
706 "boundary_correction_energy": int(energy(case["boundary_residual"]))
707 },
708 "weighted_selection_control": {
709 "coarse_detail_amplitude": 2,
710 "fine_detail_amplitude": 3,
711 "coarse_weighted_energy": 64,
712 "fine_weighted_energy": 18,
713 "raw_magnitude_selection_error": 64,
714 "weighted_selection_error": 18
715 },
716 "budget_residual_energies": [int(row[2]) for row in best_detail_budget(levels)],
717 "multiscale": {
718 "matrix_order": len(matrix),
719 "matrix_nonzero": sum(int(value != 0) for row in matrix for value in row),
720 "input_wavelet_nonzero": sum(int(value != 0) for value in input_waves),
721 "filtered_wavelet_nonzero": sum(int(value != 0) for value in output_waves),
722 "input_band_energies": [int(v) for _, v, _ in input_bands],
723 "filtered_band_energies": [int(v) for _, v, _ in output_bands],
724 "commutator_energy": int(energy(operation_commutator(x))),
725 "error_filter_coarse": int(energy(tuple(y-v for y,v in zip(case["filtered"],hp)))),
726 "error_coarse_filtered": int(energy(tuple(y-v for y,v in zip(case["filtered"],ph)))),
727 },
728 "assertions_by_group": counts,
729 "total_assertions": sum(counts.values()),
730 "scope": "Exact finite models; no infinite-domain or cross-boundary filter claim."
731 }
732 (folder/"results.json").write_text(json.dumps(record, indent=2, sort_keys=True) + "\n",
733 encoding="utf-8", newline="\n")
734
735
736 def main():
737 parser = argparse.ArgumentParser(description=__doc__)
738 parser.add_argument("--output", type=Path, help="Write reproducible CSV, SVG, and JSON records")
739 args = parser.parse_args()
740 counts = check()
741 case = default_case()
742 print("WAVELET POLYPHASE | 16 samples | circular FIR (1,2,-1)")
743 print("Original energy 416 = coarse 288 + residual 128")
744 print("Direct FIR equals both transported polyphase and Haar operators")
745 print("FIR energy 2288 = 4 * 416 + 624; two omitted transfers also cost 624")
746 print("Zero-extension boundary correction: (4,-4), squared energy 32")
747 print("Best-detail budget errors: 272, 144, 80, 48, 16, 0")
748 print("Full Haar transfer: 16 x 16 rational matrix, 56 nonzero entries")
749 print("Commutator energy 624; H(Px) error 768; P(Hx) error 880")
750 print("PASS: " + str(sum(counts.values())) + " exact assertions")
751 if args.output:
752 save(args.output, counts)
753 print("Exact tables and diagram written to " + str(args.output))
754
755
756 if __name__ == "__main__":
757 main()