1# /// script
2# dependencies = [
3#   "matplotlib",
4#   "scipy",
5# ]
6# ///
7
8"""
9Generate an FIR low-pass filter designed for an oversampled signal using the windowing method with a Kaiser window, and
10plot the frequency response.
11"""
12
13import argparse
14from math import pi
15
16import matplotlib.pyplot as plt
17import numpy as np
18from scipy.signal import firwin, freqz, kaiser_beta
19
20
21def main():
22    args = parse_args()
23
24    beta = kaiser_beta(args.sba)
25    taps = args.n * (2 * args.nz + 1)
26    coefficients = firwin(
27        taps, (args.fc / (args.cutoff_denominator_freq / 2)) / args.n, window=("kaiser", beta)
28    )
29
30    print(f"{taps} taps")
31
32    with open(args.o, "w") as f:
33        coeffs_to_write = coefficients
34        if not args.write_full:
35            coeffs_to_write = coefficients[int(len(coefficients) / 2) :]
36        for coefficient in coeffs_to_write:
37            f.write(f"{coefficient},\n")
38
39    w, h = freqz(coefficients, worN=1 << 20)
40
41    figure = plt.figure()
42    axes = figure.add_subplot(
43        title="Filter Frequency Response",
44        xlabel="Frequency (Hz)",
45        ylabel="Gain (dB)",
46        xlim=(0, min(args.fs, 4 * args.fc)),
47        ylim=(-(args.sba + 20), 20),
48    )
49    axes.grid(visible=True)
50
51    axes.plot(w * args.n * (args.fs / 2) / pi, 20 * np.log10(abs(h)))
52    axes.plot(2 * [args.fs / 2], [-300, 300])
53    axes.plot([0, args.fs], 2 * [-3])
54
55    plt.show(block=True)
56
57
58def parse_args():
59    arg_parser = argparse.ArgumentParser(
60        description="Generate a FIR low-pass filter using a Kaiser window"
61    )
62    arg_parser.add_argument(
63        "-sba",
64        default=80,
65        type=float,
66        required=False,
67        help="Stopband attenuation (dB) (default=80)",
68    )
69    arg_parser.add_argument(
70        "-n",
71        default=512,
72        type=int,
73        required=False,
74        help="Samples per zero crossing (default=512)",
75    )
76    arg_parser.add_argument(
77        "-nz", type=int, required=True, help="Number of zero crossings"
78    )
79    arg_parser.add_argument(
80        "-fs",
81        default=48000,
82        type=float,
83        required=False,
84        help="Source frequency for plotting (Hz) (default=48000)",
85    )
86    arg_parser.add_argument(
87        "-fc",
88        default=20000,
89        type=float,
90        required=False,
91        help="Cutoff frequency (Hz) (default=20000)",
92    )
93    arg_parser.add_argument(
94        "--cutoff-denominator-freq",
95        default=48000,
96        type=float,
97        required=False,
98        help="Denominator when computing cutoff value (Hz) (default=48000)"
99    )
100    arg_parser.add_argument(
101        "-o",
102        default="kaiser-fir.txt",
103        type=str,
104        required=False,
105        help="Output file (default=kaiser-fir.txt)",
106    )
107    arg_parser.add_argument(
108        "--write-full",
109        action="store_true",
110        required=False,
111        help="Write out the entire filter instead of only the right half",
112    )
113    return arg_parser.parse_args()
114
115
116if __name__ == "__main__":
117    main()