scipy_faust-iir.py

scipy_faust-iir.py

import numpy as np
from scipy.signal import iirdesign, freqz
import matplotlib.pyplot as plt

# Define the desired frequency response
def desired_frequency_response(w):
    return np.abs(np.sin(2 * np.pi * w) + 0.5 * np.cos(4 * np.pi * w))

# Frequency points
freqs = np.linspace(0, 1, 8000)
desired_response = desired_frequency_response(freqs)

# Design the filter using scipy's iirdesign
wp = 0.2  # Passband frequency
ws = 0.3  # Stopband frequency
gpass = 1  # Maximum loss in the passband (dB)
gstop = 40  # Minimum attenuation in the stopband (dB)
b, a = iirdesign(wp, ws, gpass, gstop, ftype='ellip')

# Compute the frequency response of the optimized filter
w, h = freqz(b, a, worN=8000)

# Plot the desired and actual frequency response
# plt.figure(figsize=(10, 6))
# plt.plot(freqs, desired_response, label='Desired Response', linestyle='--')
# plt.plot(w / np.pi, np.abs(h), label='Actual Response')
# plt.title('IIR Filter Design using Least-Squares Optimization')
# plt.xlabel('Normalized Frequency')
# plt.ylabel('Magnitude')
# plt.legend()
# plt.grid()
# plt.show()

print("Feedforward coefficients (b):", b)
print("Feedback coefficients (a):", a)



faust_code_template = """

import("filter.lib");

N = {};
bCoeffs = {};
aCoeffs = {};

process = _ <: iir(bCoeffs, aCoeffs);

"""



b_str = ", ".join(str(coeff) for coeff in b)
a_str = ", ".join(str(coeff) for coeff in a[1:])  # Exclude the first element of 'a'


faust_code = faust_code_template.format(b.size , b_str, a_str)


print(faust_code)

with open("optimized_iir_filter.dsp", "w") as f:
    f.write(faust_code)