Polyphase and FIR tap split / exact program
A forward-shift filter, exact signed Euler inputs, complete finite windows, and an explicit missing-advance control.
Read the worked tap split Download Python · Exact results and test ranges
Run python -B euler_polyphase_taps.py with Python 3.10 or later. Only the standard library is needed.
1 #!/usr/bin/env python3
2 """An exact forward-FIR example: split both the signal and the filter taps.
3
4 Run: python -B euler_polyphase_taps.py
5 Save the finite evidence: python -B euler_polyphase_taps.py --output results.json
6 Only complete windows of the supplied finite signal are returned. No padding,
7 periodic extension, floating arithmetic, or Fourier transform is used.
8 """
9 import argparse
10 import json
11 from math import comb
12 from pathlib import Path
13
14
15 TAPS = (1, -2, 3, 1)
16 TEST_FILTERS = (TAPS, (1,), (0, 1), (1, 0, 0, 1), (0, -2, 0, 3), (1, -1, 2, -3, 4))
17
18
19 def signed_even_euler(count):
20 values = []
21 for n in range(count):
22 values.append(1 if n == 0 else -sum(comb(2*n, 2*j)*values[j] for j in range(n)))
23 return values
24
25
26 def direct_filter(signal, taps):
27 if not taps:
28 raise ValueError('A filter must have at least one tap')
29 return [sum(h*signal[n+j] for j, h in enumerate(taps))
30 for n in range(max(0, len(signal)-len(taps)+1))]
31
32
33 def split_filter(signal, taps):
34 if not taps:
35 raise ValueError('A filter must have at least one tap')
36 a0, a1 = signal[::2], signal[1::2]
37 h0, h1 = taps[::2], taps[1::2]
38 count = max(0, len(signal)-len(taps)+1)
39 even = [sum(h*a0[q+j] for j, h in enumerate(h0))
40 + sum(h*a1[q+j] for j, h in enumerate(h1))
41 for q in range((count+1)//2)]
42 odd = [sum(h*a1[q+j] for j, h in enumerate(h0))
43 + sum(h*a0[q+j+1] for j, h in enumerate(h1))
44 for q in range(count//2)]
45 return [even[n//2] if n % 2 == 0 else odd[n//2] for n in range(count)]
46
47
48 def compute():
49 signal = signed_even_euler(16)
50 assert signal[:8] == [1, -1, 5, -61, 1385, -50521, 2702765, -199360981]
51 direct, split = direct_filter(signal, TAPS), split_filter(signal, TAPS)
52 assert direct == split
53 assert direct[:4] == [-43, 1191, -46239, 2548371]
54 even, odd = signal[::2], signal[1::2]
55 rows = []
56 for q in range(2):
57 even_part = even[q]+3*even[q+1]
58 odd_part = -2*odd[q]+odd[q+1]
59 assert even_part+odd_part == direct[2*q]
60 rows.append(dict(index=2*q, phase='even', first=even_part, second=odd_part,
61 output=direct[2*q]))
62 even_part = -2*even[q+1]+even[q+2]
63 odd_part = odd[q]+3*odd[q+1]
64 assert even_part+odd_part == direct[2*q+1]
65 rows.append(dict(index=2*q+1, phase='odd', first=even_part, second=odd_part,
66 output=direct[2*q+1]))
67 missing_advance = (-2*even[0]+even[1])+(odd[0]+3*odd[1])
68 assert missing_advance == -181 and direct[1] == 1191
69 comparisons = 0
70 for length in range(25):
71 signals = (signal[:length] if length <= 16 else signed_even_euler(length),
72 [1]*length, [(-1)**n for n in range(length)], list(range(length)))
73 for values in signals:
74 for taps in TEST_FILTERS:
75 assert split_filter(values, taps) == direct_filter(values, taps), (length, taps)
76 comparisons += 1
77 impulse_comparisons = 0
78 for length in range(1, 13):
79 for position in range(length):
80 values = [int(n == position) for n in range(length)]
81 for taps in TEST_FILTERS:
82 assert split_filter(values, taps) == direct_filter(values, taps), (length, position, taps)
83 impulse_comparisons += 1
84 for implementation in (direct_filter, split_filter):
85 try:
86 implementation([1, 2], ())
87 except ValueError:
88 pass
89 else:
90 raise AssertionError('An empty tap list was accepted')
91 return dict(status='PASS', convention='Forward shift S a(n)=a(n+1); complete finite windows only',
92 input=signal, taps=list(TAPS), even_taps=list(TAPS[::2]), odd_taps=list(TAPS[1::2]),
93 direct=direct, polyphase=split, worked_rows=rows,
94 missing_advance=dict(incorrect=missing_advance, correct=direct[1]),
95 finite_checks=dict(signal_comparisons=comparisons, impulse_comparisons=impulse_comparisons,
96 prefix_lengths=[0, 24], impulse_lengths=[1, 12],
97 filters=[list(t) for t in TEST_FILTERS], empty_filter_rejected=True),
98 scope='Exact finite comparisons and the displayed arithmetic. The general phase law is derived in the worked page; these tests are not its proof.')
99
100
101 def main():
102 parser = argparse.ArgumentParser(description=__doc__)
103 parser.add_argument('--output', type=Path)
104 args = parser.parse_args()
105 report = compute()
106 if args.output:
107 if args.output.exists():
108 raise FileExistsError('Choose a new result file')
109 args.output.parent.mkdir(parents=True, exist_ok=True)
110 args.output.write_text(json.dumps(report, indent=2)+'\n', encoding='utf-8', newline='\n')
111 print('Even taps: 1, 3')
112 print('Odd taps: -2, 1')
113 for label, key in (('Direct outputs', 'direct'), ('Polyphase outputs', 'polyphase')):
114 print(label+': '+', '.join(map(str, report[key][:4])))
115 print('Missing-advance control: -181 (correct: 1191)')
116 print('Finite comparisons: 600 signal cases and 468 impulse cases')
117 print('PASS: exact polyphase/tap split')
118
119
120 if __name__ == '__main__':
121 main()