Bode plot of a 157 order transfer function
조회 수: 1 (최근 30일)
이전 댓글 표시
I am trying to evaluate a cascade of 480 stages of second-order low-pass filters (to mimic human cochleas) in s-domain. The resulted filter is the 157 order (has 158 non-zero coeffiecients in denominator). The problem is that I got nothing from functions bode(myfilter) or lsim(myfilter, input_signal, t). But I can plot some kind of frequency response with freqresp(). Is it because bode() and lsim() have order limitations? I also attached the code below:
% A second order low pass filter
% Cascade of 480 stages
Q = 0.79;
numerator = 1;
mResult = 1; % Accumulated transfer function
for k = 1:480
tau(k) = 1.014^(k-1)/(40000*pi);
myFilters(k) = tf(1, [tau(k)^2, 1/Q*tau(k), 1]);
end
% Place selection
for t = 1:480
mResult = mResult*myFilters(t);
end
f = 1:1:30000; % 1-30kHz
w = 2*pi*f;
filterUndertest = mResult;
out = freqresp(filterUndertest, w);
for k = 1:30000
AbsOut (k) = abs(out(:, :, k));
end
semilogx(f,AbsOut); % This thing works
bode(mResult); % But this not
댓글 수: 0
채택된 답변
Paul
2023년 10월 25일
편집: Paul
2023년 10월 26일
Hi Dongxu Guo,
That's quite a filter!
Using a tf for a for such a high order filter is not likely to work due to large rounding errors. BTW, isn't the filter order 480*2 = 960?
For Bode plots, my quick experiments indicate the frd representation seems to be the best approach, though I have no idea if the final result looks like it should.
% A second order low pass filter
% Cascade of 480 stages
Q = 0.79;
numerator = 1;
f = 1:1:30000; % 1-30kHz
f = logspace(0,5,1000); % use less points for run time
f(f>30000) = [];
w = 2*pi*f;
mResult = frd(tf(1),w); % Accumulated transfer function
N = 480;
for k = 1:N %1:480
tau(k) = 1.014^(k-1)/(40000*pi);
myFilters{k} = tf(1, [tau(k)^2, 1/Q*tau(k), 1]);
end
% Place selection
for t = 1:N %480
mResult = mResult*frd(myFilters{t},w);
end
%{
f = 1:1:30000; % 1-30kHz
w = 2*pi*f;
filterUndertest = mResult;
out = freqresp(filterUndertest, w);
for k = 1:30000
AbsOut (k) = abs(out(:, :, k));
end
%}
[m,p,w]=bode(mResult,w);
bode(mResult)
The Bode plot doesn't go all the way to f(end) because eventually the magnitudes are numerically 0 (I assume the phase plot stops to be consistent with the magnitude plot, though the calculated phase wil likely be inncacurate anyway).
m(end),p(end)
figure
semilogx(f,abs(m(:)));
For lsim, you might want to try using a loop to lsim each one with the output of the previous stage being the input of the current stage in the loop.
댓글 수: 3
Paul
2023년 10월 26일
Here is the code from the question
Q = 0.79;
numerator = 1;
mResult = 1; % Accumulated transfer function
for k = 1:480
tau(k) = 1.014^(k-1)/(40000*pi);
myFilters(k) = tf(1, [tau(k)^2, 1/Q*tau(k), 1]);
end
I just want to point out that myFilters constructed this way is not really an arrary of transfer functions, but rather a transfer function with 1 output and 480 inputs. Assuming we really want to just store separate transfer functions, consider using a cell array (thought the current code seems to work fine), which I've now updated in the Answer.
whos myFilters
size(myFilters)
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Digital Filter Analysis에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!