Bode plot of a 157 order transfer function
1 view (last 30 days)
Show older comments
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 Comments
Accepted Answer
Paul
on 25 Oct 2023
Edited: Paul
on 26 Oct 2023
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 Comments
Paul
on 26 Oct 2023
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)
More Answers (0)
See Also
Categories
Find more on Digital Filter Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!