Fast implementation of calculating peaks/local maxima of magnitude of transfer function

Hi,

I need to calculate different peak values of the magnitude of the transfer function. The peaks are all local maxima, not just the highest value. I have the code below to illustrate an example. You could see that the values in pks_vect_maxima are matching with the values in Bode plot. This is my current implementation and it is slow for my application. I think going through Bode route makes it slow. Is there a better/faster way for implementation? I am looking for a solution that works for general cases of transfer function order.

Thanks,

clear all;
options_bode = bodeoptions;
options_bode.MagUnits = 'abs';
G_tf = tf([0.2461,   -0.0829,   -0.0673,    0.0829,   -0.1788], ...
    [1.0000,   -1.2841,    1.0272,   -0.9329,    0.5279],1);
figure;
hold on;
bodemag(G_tf,'g',options_bode)
% Calculate local maxima/peaks value of transfer function
[mag,phase,wout] = bode(G_tf);
mag_vect = mag(:);
[pks_vect_maxima,locs_maxima] = findpeaks(mag_vect);

Answer 1

[Star Strider]

Using the Signal Processing Toolbox function is significantly faster:

b = [0.2461,   -0.0829,   -0.0673,    0.0829,   -0.1788];
a = [1.0000,   -1.2841,    1.0272,   -0.9329,    0.5279];

figure
[h,w] = freqz(b, a, 2^14);
[pks_vect_maxima2,locs_maxima3] = findpeaks(abs(h));

with the bode and findpeaks calls requiring about 1.7 seconds and the freqz and findpeaks calls requiring about 0.05 seconds when I tested them

You need to decide if the extrra coding required is worth the effort.

My test code:

options_bode = bodeoptions;
options_bode.MagUnits = 'abs';
G_filt_discrete = tf([0.2461,   -0.0829,   -0.0673,    0.0829,   -0.1788], ...
    [1.0000,   -1.2841,    1.0272,   -0.9329,    0.5279],1);
t1 = clock;
figure;
hold on;
bodemag(G_filt_discrete,'g',options_bode)
% Calculate local maxima/peaks value of transfer function
[mag,phase,wout] = bode(G_filt_discrete);
mag_vect = mag(:);
[pks_vect_maxima,locs_maxima] = findpeaks(mag_vect);
t2 = clock;
pkw = wout(locs_maxima)

b = [0.2461,   -0.0829,   -0.0673,    0.0829,   -0.1788];
a = [1.0000,   -1.2841,    1.0272,   -0.9329,    0.5279];

figure
[h,w] = freqz(b, a, 2^14);
[pks_vect_maxima2,locs_maxima2] = findpeaks(abs(h));
t3 = clock
pkw2 = w(locs_maxima2)
et1 = etime(t2,t1)
et2 = etime(t3,t2)

[Dave Lee]

Hi Star Strider: thanks for the input

It is interesting that if I run one time (for-loop in the code below runs one iteration), freqz (0.03sec) is faster than bode (0.18sec). When for-loop iterations are 10000times, freqz (27.3sec) becomes slower than bode (24.5sec). For my main problem, how slower freqz than bode is even more significant.

clear all;
clc;
b = [0.2461,   -0.0829,   -0.0673,    0.0829,   -0.1788];
a = [1.0000,   -1.2841,    1.0272,   -0.9329,    0.5279];
G_filt_discrete = tf(b,a,1);
tic
for xx=1:10000
    [mag,phase,wout] = bode(G_filt_discrete);
    mag_vect = mag(:);
    [pks_vect_maxima,locs_maxima] = findpeaks(mag_vect);
end
toc
tic;
for xx=1:10000
    [h,w] = freqz(b, a, 2^14);
    [pks_vect_maxima2,locs_maxima2] = findpeaks(abs(h));
end
toc

[Dave Lee]

I tried reducing the length and it helped for making it faster than bode for the for-loop case above. But it is not for another case.

I think comparison between two is application-specific. And this alternative is helpful for others to see an option to consider. If possible, please leave your answer here.

Thanks,

[Star Strider]

That is the only alternative I can suggest. You can experiment with decreasing the fft lengths (from 2^14 to 2^10 or something shorter). I doubt there any faster possibilities.

The tradeoff is speed or precision in calculating the function from which to derive the peak amplitudes and locations. Longer fft vectors are more precise (providing increased frequency resolution), and are by definition slower. It depends on what you want.

I will delete my Answer in a few hours, since it did not turn out to be helpful.