-
Notifications
You must be signed in to change notification settings - Fork 0
/
fir_filter_ls.py
134 lines (100 loc) · 3.49 KB
/
fir_filter_ls.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Synthesise FIR filter from complex valued frequency response samples.
@author: Arnfinn Eielsen
@date: 08.04.2024
@license: BSD 3-Clause
"""
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy import signal
import control as ct
def fir_filter_ls(H, p):
"""
Synthesise FIR filter from complex valued frequency response samples,
using least-squares. No restriction on phase response.
Inputs:
H - frequency response samples from 0 to 2*pi
p - desired FIR filter length
Outputs:
alpha - FIR filter coefficients
alpha_win - windows FIR filter coefficients
beta - denominator that can be used to null the phase when computing the frequency response
"""
M = H.size # no. of frequency samples
X = np.zeros((M,p), dtype=np.complex_)
if np.mod(p,2): # odd
q = int((p-1)/2)
else: # even
q = int(p/2)
K = 2*np.pi/M
ps = np.arange(0,p) #0:p-1
pq = K*(ps-q)
for k in range(0, M): # iterate over frequency samples
X[k,:] = np.exp(-1j*k*pq)
b = H.reshape(-1, 1) # column vector
# min(||b - X*alpha||)
alpha = np.linalg.lstsq(X, b, rcond=None)[0] # FIR coefficients by least-squares
alpha = np.real(alpha.reshape(1, -1)) # real valued row vector
win = np.hanning(alpha.size)
alpha_win = win*alpha # reduce ripples due to implicit rectangular windowing
beta = np.r_[1.0, np.zeros(q-1)]
return alpha.squeeze(), alpha_win.squeeze(), beta.squeeze()
def main():
"""
Test the method.
"""
M = int(1e4) # no. of frequency samples
w = np.linspace(0.0, 2*np.pi, M) # sample whole circle
match 1:
case 1:
sys = ct.drss(5, outputs=1, inputs=1) # random, stable LTI system
Hss = signal.dlti(sys.A, sys.B, sys.C, sys.D, dt=1)
wH_fr, H_fr = signal.dfreqresp(Hss, w)
case 2:
hu = np.linspace(0, 1, int(w.size/2))
hd = np.linspace(1, 0, int(w.size/2))
H_fr = np.r_[hu, hd]
case 3:
hs = 0.0125*np.ones(int(1000))
hp = 3*np.ones(int(8000))
H_fr = np.r_[hs, hp, hs]
I = np.argwhere(w < np.pi)
plt.plot(w, abs(H_fr))
plt.show()
p = 200 # filter length
alpha, alpha_win, beta = fir_filter_ls(H_fr, p)
wh, h = signal.freqz(alpha_win.squeeze(), beta, w)
#fig, ax1 = plt.subplots()
#ax1.set_title('Digital filter frequency response')
# ax1.plot(w, 20*np.log10(h), 'b')
#ax1.plot(wh, abs(h), 'b')
plt.plot(wh, abs(h), 'b')
# ax1.set_ylabel('Amplitude [dB]', color='b')
# ax1.set_xlabel('Frequency [rad/sample]')
# ax2 = ax1.twinx()
# angles = np.unwrap(np.angle(h))
# ax2.plot(wh, angles, 'g')
# ax2.set_ylabel('Angle (radians)', color='g')
# ax2.grid(True)
# ax2.axis('tight')
plt.show()
"""
switch nargin
case 2
case 1
p = 50; % filter length
otherwise
HM = drss(5); % random LTI system prototype/reference
HMtf = tf(HM);
bref = HMtf.Numerator{:};
aref = HMtf.Denominator{:};
M = 1e3; % no. of frequency samples
w = linspace(0,2*pi,M); % sample whole circle
H = freqz(bref,aref,w); % frequency response samples
p = 50; % filter length
end
"""
if __name__ == "__main__":
main()