Converting filter equation to frequency response
Note: Part of the problem is solved but I've kept the original post for clarity.
I was told this "should be easy" for a filter designer. I am not so sure though, but I'll ask anyway :)
What I need to construct from these equations is a function that takes a frequency range and produces the integrated response of the filter across that range:
float filter_response(float f_low, float f_high) {
return {magic};
}
The frequency response of the filter is known:
And the expected output from desired function is known from a table for verification:
| š¯‘“c | š¯‘“low | š¯‘“high | output |
|---|---|---|---|
| 4 | 3.5 | 4.5 | 0.375 |
| 5 | 4.4 | 5.6 | 0.545 |
| 6.3 | 5.6 | 7.0 | 0.727 |
| 8 | 7.1 | 8.9 | 0.873 |
| 10 | 8.8 | 11.2 | 0.951 |
| 12.5 | 11.1 | 13.9 | 0.958 |
| 16 | 14.2 | 17.8 | 0.896 |
| 20 | 17.7 | 22.3 | 0.782 |
| 25 | 22.1 | 27.9 | 0.647 |
| 31.5 | 27.9 | 35.1 | 0.519 |
| 40 | 35.4 | 44.6 | 0.411 |
| 50 | 44.2 | 55.8 | 0.324 |
| 63 | 55.8 | 70.2 | 0.256 |
| 80 | 70.8 | 89.2 | 0.202 |
| 100 | 88.5 | 111.5 | 0.16 |
| 125 | 110.6 | 139.4 | 0.127 |
| 160 | 141.6 | 178.4 | 0.101 |
| 200 | 177.0 | 223.0 | 0.0799 |
| 250 | 221.2 | 278.8 | 0.0634 |
| 315 | 278.8 | 351.2 | 0.0503 |
| 400 | 354.0 | 446.0 | 0.0398 |
| 500 | 442.5 | 557.5 | 0.0315 |
| 630 | 557.5 | 702.5 | 0.0245 |
| 800 | 708.0 | 892.0 | 0.0186 |
| 1000 | 885.0 | 1115.0 | 0.0135 |
| 1250 | 1106.2 | 1393.8 | 0.00894 |
| 1600 | 1416.0 | 1784.0 | 0.00538 |
| 2000 | 1770.0 | 2230.0 | 0.00295 |
Where š¯‘“c is the center frequency, š¯‘“low-š¯‘“high the frequency range, and output the expected output. Only š¯‘“c was originally transcribed from the table: š¯‘“low-š¯‘“high were calculated by me and represent the frequency range of one third octave for the given "mid"-frequency. I provided these for clarity as they directly relate to the desired function.
This is the data that I have and while I am sure this can be done I also understand that this might be beyond what could be expected to get help with online (unless it really is as easy as it was presented to me).
Regardless, any pointers in the right direction would be much appreciated.
Thank you <3
Since I am neither a filter designer or very well versed in math, I've mainly looked at extracting data from the graph or by interpolating the given values in the table. While this works for now it does not give me the accuracy or trustworthiness I am looking for. I read code in pretty much any language and I understand FIR and IIR filters fairly well.
I've tried designing a filter with Python, SciPy, Numpy based on the given filter constants to see if I could replicate the given response curve, but I've been unsuccessful.
I've also tried AI resources to help me cut through these formulas to give me something that describes the filter in terms that fit my limited knowledge in the field, but they appear wholly unable to deal with this problem.
In the end, I expect to be able to solve this in several different ways. While a function numerically calculating this would be fantastic I could, in theory, also analyze an FFT/FIR filter to produce a new table with the range and precision needed.
Update:
import math
import numpy as np
import matplotlib.pyplot as plt
# constants
pi, f_1, f_2, Q_1, f_3, f_4, Q_2, K = math.pi, 6.31, 1258.9, 0.71, 15.915, 15.915, 0.64, 1
# band-limiting transfer function
def band_limit(s):
return ((s**2)*4*(pi**2)*(f_2**2))/(((s**2)+2*pi*f_1*s/Q_1+4*(pi**2)*(f_1**2))*((s**2)+2*pi*f_2*s/Q_1+4*(pi**2)*(f_2**2)))
# frequency-weighting transfer function
def frequency_weight(s):
return ((s+2*pi*f_3)*2*pi*K*(f_4**2))/(((s**2)+2*pi*f_4*s/Q_2+4*(pi**2)*(f_4**2))*f_3)
px, py, py1, py2 = [], [], [], []
for x in np.arange(1,1000,0.1):
w = x*2*pi # angular frequency
px.append(x)
py1.append(band_limit(w))
py2.append(frequency_weight(w))
py.append(band_limit(w)*frequency_weight(w))
plt.plot(px, py1, label = "band-limit filter")
plt.plot(px, py2, label = "frequency-weighting filter")
plt.plot(px, py, label = "total")
plt.yscale("log");
plt.xscale("log");
plt.ylim(0.01, 1);
plt.xlim(1, 1000);
plt.legend()
plt.show()
This python program produces the following plot:
I've overlaid the original plot for comparison.
The plots do not match, but they appear not entirely dissimilar which tells me I am on the right path. Is it possible to use a real-only valued function to calculate the amplitude response at each frequency or is the error I am getting due to not performing the calculations with complex numbers?
Update:
Apparently the calculation must be done with complex math as the following python code and diagram shows (a near perfect match):
import math
import numpy as np
import matplotlib.pyplot as plt
import cmath
# constants
pi, f_1, f_2, Q_1, f_3, f_4, Q_2, K = math.pi, 6.31, 1258.9, 0.71, 15.915, 15.915, 0.64, 1
# band-limiting transfer function
def band_limit(s):
z = ((s**2)*4*(pi**2)*(f_2**2))/(((s**2)+2*pi*f_1*s/Q_1+4*(pi**2)*(f_1**2))*((s**2)+2*pi*f_2*s/Q_1+4*(pi**2)*(f_2**2)))
return math.sqrt(z.real**2+z.imag**2)
# frequency-weighting transfer function
def frequency_weight(s):
z = ((s+2*pi*f_3)*2*pi*K*(f_4**2))/(((s**2)+2*pi*f_4*s/Q_2+4*(pi**2)*(f_4**2))*f_3)
return math.sqrt(z.real**2+z.imag**2)
px, py, py1, py2 = [], [], [], []
for x in np.arange(1,1259.9,0.1):
w = x*2*pi # angular frequency
px.append(x)
py1.append(band_limit(complex(1.0,w)))
py2.append(frequency_weight(complex(1.0,w)))
py.append(band_limit(complex(1.0,w))*frequency_weight(complex(1.0,w)))
plt.plot(px, py1, label = "band-limit filter")
plt.plot(px, py2, label = "frequency-weighting filter")
plt.plot(px, py, label = "total")
plt.yscale("log");
plt.xscale("log");
plt.ylim(0.01, 1);
plt.xlim(1, 1259.9);
plt.legend()
plt.show()
The only question left is how do I integrate this across a range to produce the function originally asked for, such that it could produce the verification table?
Seems fine, though I don't know what your tolerances are:
from io import StringIO
import numpy as np
import pandas as pd
import scipy.integrate
from matplotlib import pyplot as plt
def Hb(
s: float | np.ndarray,
f1: float = 6.310,
f2: float = 1258.9,
Q1: float = 0.71,
) -> np.ndarray:
from numpy import pi
w1 = 2*pi*f1
w2 = 2*pi*f2
return (s * w2)**2 / (
s**2 + w1*s/Q1 + w1**2
) / (
s**2 + w2*s/Q1 + w2**2
)
def Hw(
s: float | np.ndarray,
f3: float = 15.915,
f4: float = 15.915,
Q2: float = 0.64,
K: int = 1,
) -> np.ndarray:
from numpy import pi
w3 = 2*pi*f3
w4 = 2*pi*f4
return (s + w3) * (
K * w4**2
) / (
s**2 + w4*s/Q2 + w4**2
) / w3
def H(f: float | np.ndarray) -> np.ndarray:
s = 2j * np.pi * f
return Hb(s) * Hw(s)
def load_reference() -> pd.DataFrame:
content = '''
fc flow fhigh output
4 3.5 4.5 0.375
5 4.4 5.6 0.545
6.3 5.6 7.0 0.727
8 7.1 8.9 0.873
10 8.8 11.2 0.951
12.5 11.1 13.9 0.958
16 14.2 17.8 0.896
20 17.7 22.3 0.782
25 22.1 27.9 0.647
31.5 27.9 35.1 0.519
40 35.4 44.6 0.411
50 44.2 55.8 0.324
63 55.8 70.2 0.256
80 70.8 89.2 0.202
100 88.5 111.5 0.16
125 110.6 139.4 0.127
160 141.6 178.4 0.101
200 177.0 223.0 0.0799
250 221.2 278.8 0.0634
315 278.8 351.2 0.0503
400 354.0 446.0 0.0398
500 442.5 557.5 0.0315
630 557.5 702.5 0.0245
800 708.0 892.0 0.0186
1000 885.0 1115.0 0.0135
1250 1106.2 1393.8 0.00894
1600 1416.0 1784.0 0.00538
2000 1770.0 2230.0 0.00295
'''
with StringIO(content) as f:
return pd.read_csv(f, delim_whitespace=True)
def get_actual(ref: pd.DataFrame) -> np.ndarray:
return np.abs(H(ref['fc']))
def get_interval_mean(ref: pd.DataFrame) -> np.ndarray:
# This is needed because scipy can't vectorise bounds
integ = np.array([
scipy.integrate.quad(
func=H, a=row['flow'], b=row['fhigh'], complex_func=True,
)[0]
for idx, row in ref.iterrows()
])
means = integ/(ref['fhigh'] - ref['flow'])
return means.abs()
def verify_band_limiting(ref: pd.DataFrame) -> None:
assert np.allclose(
a=ref['actual'],
b=ref['output'],
rtol=0.05, atol=0,
)
def plot(ref: pd.DataFrame) -> plt.Figure:
fig, (ax_transfer, ax_err) = plt.subplots(ncols=2)
f_actual = np.exp(np.linspace(np.log(1), np.log(2000), 250))
ax_transfer.loglog(f_actual, np.abs(H(f_actual)), label='actual', color='orange')
ax_transfer.scatter(
ref['fc'], ref['output'],
label='reference', color='black', s=1, zorder=10,
)
ax_transfer.set_xlabel('Frequency (Hz)')
ax_transfer.set_ylabel('Magnitude')
legend_kwargs = {'label': 'interval mean'}
for idx, row in ref.iterrows():
ax_transfer.plot(
[row['flow'], row['fc'], row['fhigh']],
[row['interval_mean'], row['actual'], row['interval_mean']],
color='lightgreen', **legend_kwargs,
)
legend_kwargs = {}
ax_transfer.legend()
ax_err.semilogx(
ref['fc'], np.zeros_like(ref['actual']),
label='actual', color='orange')
ax_err.scatter(
ref['fc'], ref['output']/ref['actual'] - 1,
label='reference', color='black', marker='+')
ax_err.scatter(
ref['fc'], ref['interval_mean']/ref['actual'] - 1,
label='interval_mean', color='green', marker='+')
ax_err.set_xlabel('Frequency (Hz)')
ax_err.set_ylabel('Relative magnitude error')
ax_err.legend()
return fig
def main() -> None:
ref = load_reference()
ref['actual'] = get_actual(ref)
ref['interval_mean'] = get_interval_mean(ref)
verify_band_limiting(ref)
print(ref)
plot(ref=ref)
plt.show()
if __name__ == '__main__':
main()
fc flow fhigh output actual interval_mean
0 4.0 3.5 4.5 0.37500 0.379550 0.378397
1 5.0 4.4 5.6 0.54500 0.545093 0.540838
2 6.3 5.6 7.0 0.72700 0.728996 0.721519
3 8.0 7.1 8.9 0.87300 0.879784 0.870259
4 10.0 8.8 11.2 0.95100 0.954051 0.943934
5 12.5 11.1 13.9 0.95800 0.960305 0.952859
6 16.0 14.2 17.8 0.89600 0.893056 0.887968
7 20.0 17.7 22.3 0.78200 0.781253 0.779000
8 25.0 22.1 27.9 0.64700 0.650180 0.650352
9 31.5 27.9 35.1 0.51900 0.521334 0.522628
10 40.0 35.4 44.6 0.41100 0.409186 0.410701
11 50.0 44.2 55.8 0.32400 0.325219 0.326599
12 63.0 55.8 70.2 0.25600 0.256466 0.257566
13 80.0 70.8 89.2 0.20200 0.200944 0.201825
14 100.0 88.5 111.5 0.16000 0.160218 0.160919
15 125.0 110.6 139.4 0.12700 0.127881 0.128439
16 160.0 141.6 178.4 0.10100 0.099738 0.100168
17 200.0 177.0 223.0 0.07990 0.079705 0.080044
18 250.0 221.2 278.8 0.06340 0.063703 0.063969
19 315.0 278.8 351.2 0.05030 0.050487 0.050690
20 400.0 354.0 446.0 0.03980 0.039636 0.039786
21 500.0 442.5 557.5 0.03150 0.031490 0.031595
22 630.0 557.5 702.5 0.02450 0.024557 0.024621
23 800.0 708.0 892.0 0.01860 0.018501 0.018531
24 1000.0 885.0 1115.0 0.01350 0.013510 0.013529
25 1250.0 1106.2 1393.8 0.00894 0.009104 0.009141
26 1600.0 1416.0 1784.0 0.00538 0.005255 0.005314
27 2000.0 1770.0 2230.0 0.00295 0.002939 0.002992
- Math.Sin() gives incorrect value
- How to run my python script when the sunOS is start booting
- Express-session: not resetting cookie expiration on each request
- Getting a stack overflow exception when normalizing a vector
- Edit default summary function in R gives error for multiple variables
- What was a For loop? Why isn't it needed in R?
- How to use download button in shiny and save results in various formats (csv, texte, pdf, spss...)?
- Why are there two assignment operators, `<-` and `->` in R?
- lm()$assign: what is it?
- How to get the value of list(...) in R and S functions
- Design matrix for MLM from library(lme4) with fixed and random effects
- how to generate elements not included in my sample
- Create a matrix with gradually changing values without a for loop
- Emacs ESS and S-plus ( S+ ) 8.1 compatability
- How to lag date-index in a time-series in R?
- Nonlinear regression in R / S
- Calling R from S-Plus?