Wiemy, że ogólnie funkcję transferu filtra zapewniają:
H(z)=∑Mk=0bkz−k∑Nk=0akz−k
Teraz zamień z=ejω aby ocenić funkcję przenoszenia na okręgu jednostkowym:
H(ejω)=∑Mk=0bke−jωk∑Nk=0ake−jωk
Staje się to zatem tylko problemem oceny wielomianowej przy danym ω . Oto kroki:
- Utwórz wektor częstotliwości kątowych ω=[0,…,π] dla pierwszej połowy widma (nie trzeba zwiększać do 2π ) i zapisz go
w
.
- Przed obliczeniem wykładnika e−jω na wszystkich z nich i zapisz go w zmiennej
ze
.
- Użyj
polyval
funkcji, aby obliczyć wartości licznika i mianownika, wywołując polyval(b, ze)
:, podziel je i zapisz H
. Ponieważ interesuje nas amplituda, weź bezwzględną wartość wyniku.
- Przelicz na skalę dB za pomocą: HdB=20log10H - w tym przypadku 1 jest wartością odniesienia.
Umieszczenie tego wszystkiego w kodzie:
%% Filter definition
a = [1 -0.5 -0.25]; % Some filter with lot's of static gain
b = [1 3 2];
%% My freqz calculation
N = 1024; % Number of points to evaluate at
upp = pi; % Evaluate only up to fs/2
% Create the vector of angular frequencies at one more point.
% After that remove the last element (Nyquist frequency)
w = linspace(0, pi, N+1);
w(end) = [];
ze = exp(-1j*w); % Pre-compute exponent
H = polyval(b, ze)./polyval(a, ze); % Evaluate transfer function and take the amplitude
Ha = abs(H);
Hdb = 20*log10(Ha); % Convert to dB scale
wn = w/pi;
% Plot and set axis limits
xlim = ([0 1]);
plot(wn, Hdb)
grid on
%% MATLAB freqz
figure
freqz(b,a)
Oryginalna produkcja freqz
:
Wyjście mojego skryptu:
I szybkie porównanie w skali liniowej - wygląda świetnie!
[h_f, w_f] = freqz(b,a);
figure
xlim = ([0 1]);
plot(w, Ha) % mine
grid on
hold on
plot(w_f, abs(h_f), '--r') % MATLAB
legend({'my freqz','MATLAB freqz'})
Teraz możesz przepisać go na jakąś funkcję i dodać kilka warunków, aby była bardziej użyteczna.
Innym sposobem (wcześniej zaproponowanym jest bardziej wiarygodny) byłoby wykorzystanie podstawowej właściwości, że odpowiedź częstotliwościowa filtra jest transformatą Fouriera jego odpowiedzi impulsowej:
H(ω)=F{h(t)}
Dlatego musisz wprowadzić sygnał do swojego systemu δ(t) , obliczyć odpowiedź filtra i wziąć FFT:
d = [zeros(1,length(w_f)) 1 zeros(1,length(w_f)-1)];
h = filter(b, a, d);
HH = abs(fft(h));
HH = HH(1:length(w_f));
Dla porównania spowoduje to:
a = [1 -0.5 -0.25]; % Some filter with lot's of static gain
. Czy możesz mi wyjaśnić tutaj wybór tych parametrów. Czytam instrukcję mojego Matlaba i jest napisana[h,w] = freqz(hfilt,n)
w jednej części synapsy. Dajesz dwa filtry (a, b) dofreqz
. Czy oba filtry są włączonehfilt
? Lub jeden wn
?jest to tylko dodatek do odpowiedzi Jojeka, który jest bardziej ogólny i doskonale dobry, gdy używana jest matematyka o podwójnej precyzji. gdy występuje mniejsza precyzja, pojawia się „problem kosinusowy”, który pojawia się, gdy albo częstotliwość w odpowiedzi częstotliwościowej jest bardzo niska (znacznie niższa niż Nyquist), a także gdy częstotliwości rezonansowe filtra są bardzo niskie.
rozważ tę tożsamość trig:
który ma złożoną charakterystykę częstotliwościową
który ma kwadratową wielkość:
używając powyższej tożsamości trig, otrzymujesz do kwadratu wielkości:
whereϕ≜sin2(ω2)
if your gear is intending to plot this as dB, it comes out as
so your division turns into subtraction, but you have to be able to compute logarithms to some base or another. numerically, you will have much less trouble with this for low frequencies than doing it the apparent way.
źródło