Determine -3db gain on Freq- gain plot

10 ビュー (過去 30 日間)
Moussa Ihab
Moussa Ihab 2022 年 11 月 18 日
コメント済み: Star Strider 2022 年 11 月 18 日
Hello,
I got some freq-gain data from a filter and I am trying to determine -3db gain on the plot and get its frequency.
I am trying "interp1" but it comes with only one point.
How can I get the other one ?
Thank you
freq=t1.f;
noise=t1.noise;
gain=t1.gain;
figure
plot(freq,noise)
xlabel('Frequency')
ylabel('noise')
figure
[Max_gain,index1]=max(gain);
freq_max=freq(index1);
Gain_3db=Max_gain-3;
freq_3db=interp1(gain,freq,Gain_3db)
plot(freq,gain,'b')
ylabel('gain')
xlabel('Frequency')
hold on
plot(freq_3db,Gain_3db,'r*')

採用された回答

Star Strider
Star Strider 2022 年 11 月 18 日
編集済み: Star Strider 2022 年 11 月 18 日
Use:
find(diff(sign(gain-Gain_3db)))
or something similar, to find the indices of all the points where the curve equals the specified value, then include a short range (I use [-1 0 1]) so that interp1 can interpolate over a limited set of values. This prevents problems with functions that can vary significantly, especially those that change frequently.
This approach works regardless of the number of times the transfer function equals the desired decibel value (so it would work for a filter with multiple passbands and stopbands), so long that at least one does.
I do not have your data, however this should work as an illustration —
freq = linspace(0, 10, 150).'; % Assuming Column Vectors
gain = mag2db(exp(-(freq-5).^2)); % Create Function, Copnvert To 'dB'
[Max_gain,index1]=max(gain);
freq_max=freq(index1);
Gain_3db=Max_gain-3;
idxr = find(diff(sign(gain-Gain_3db)))+(-1:1) % Matrix Of Index Ranges
idxr = 2×3
65 66 67 83 84 85
for k = 1:size(idxr,1)
freq_3db(k) = interp1(gain(idxr(k,:)),freq(idxr(k,:)),Gain_3db);
end
freq_3db % Desired Frequencies
freq_3db = 1×2
4.4121 5.5879
figure
hp1 = plot(freq, gain, 'DisplayName','Transfer Function');
hold on
hp2 = plot(freq_3db, Gain_3db, '+r', 'DisplayName','Transfer Function -3 dB Points');
hold off
grid
legend([hp1; hp2(1)], 'Location','best')
ylim([-20, 5])
Make appropriate changes to work with your data.
EDIT — (18 Nov 2022 at 13:15)
Clarified description. Code unchanged.
EDIT — (18 Nov 2022 at 19:48)
With the provided data —
freq = [2400 2470 2540 2610 2680 2750 2820 2890 2960 3030 3100 3170 3240 3310 3380 3450 3520 3590 3660 3730 3800 3870 3940 4010 4080 4150 4220 4290 4360 4430 4500 4570 4640 4710 4780 4850 4920 4990 5060 5130 5200];
gain = [-18.6680 -14.4060 -10.2174 -6.4374 -3.7861 -2.4066 -1.7975 -1.6368 -1.7193 -1.7336 -1.7873 -1.4274 -1.7189 -1.3573 -1.6372 -1.3042 -1.2758 -1.1481 -1.0076 -1.1632 -1.0454 -1.1071 -1.1686 -1.1704 -1.3548 -1.2611 -1.4895 -1.3286 -1.6556 -1.4886 -2.1840 -2.1738 -3.9351 -3.9649 -6.5240 -7.0555 -10.1172 -11.1532 -13.9887 -16.3625 -18.6674];
freq = freq(:);
gain = gain(:);
[Max_gain,index1]=max(gain);
freq_max=freq(index1);
Gain_3db=Max_gain-3;
idxr = find(diff(sign(gain-Gain_3db)))+(-1:1) % Matrix Of Index Ranges
idxr = 2×3
3 4 5 33 34 35
for k = 1:size(idxr,1)
freq_3db(k) = interp1(gain(idxr(k,:)),freq(idxr(k,:)),Gain_3db);
end
fprintf('-3 dB Frequency = %.3f Hz\n',freq_3db)
-3 dB Frequency = 2674.152 Hz -3 dB Frequency = 4711.168 Hz
% freq_3db % Desired Frequencies
figure
hp1 = plot(freq, gain, 'DisplayName','Transfer Function');
hold on
hp2 = plot(freq_3db, Gain_3db, '+r', 'DisplayName','Transfer Function -3 dB Points');
hold off
grid
legend([hp1; hp2(1)], 'Location','best')
ylim([-20, 5])
txtc = compose(' \\leftarrow %9.3f Hz',freq_3db);
text(freq_3db+15, [1 1]*Gain_3db, txtc, 'Horiz','left', 'Vert','middle', 'Rotation',-90)
Make appropriate changes to get different results.
.
  2 件のコメント
Star Strider
Star Strider 2022 年 11 月 18 日
As always, my pleasure!

サインインしてコメントする。

その他の回答 (4 件)

cr
cr 2022 年 11 月 18 日
interp1() outputs linearly interpolated values that you are asking for. if Gain_3db is a scalar value the output is a scalar. The output you get this way gives you the freq at which gain fell by 3dB from its max value.
Perhaps the other value you are looking for is the start freq.
This is also not a great way to do what you seem to be looking for. It only works if gain values are monotonic.

Askic V
Askic V 2022 年 11 月 18 日
I think this code can help you, it is something I often use. It is not most efficient, but it works.
Please have a look at the example:
clear
clc
close all
num = 1;
den = [1 1];
g = tf(num,den);
[mag, phase, w] = bode(g);
magdb = squeeze(20*log10(mag));
semilogx(w,magdb);
ylim([-15, 0]);
grid on;
ind = find(magdb>=-3,1,'last');
% find two nearest points
P1w = [w(ind), w(ind+1)];
P2m = [magdb(ind), magdb(ind+1)];
% Find slope and intercept
c = [[1; 1] P1w(:)]\P2m(:);
% y = ax+b;
a = c(2);
b = c(1);
% Find angular frequency where magdb is -3 dB
w_3db = (-3-b)/a;
% Check if this is the case
mag_3db = a*w_3db+b;
hold on
plot(w_3db, mag_3db,'ro');
text(w_3db, mag_3db, [' -3 dB point,', newline, 'w = ', num2str(w_3db), 'rad/s'])

Moussa Ihab
Moussa Ihab 2022 年 11 月 18 日
Thanks for all answers you wrote.
I might have a small data which afftects plot and interceptions.
I have tried
freq_H=interp1(gain,freq,Gain_3db,)
freq_L=interp1(gain,freq,Gain_3db, 'previous')
and I got this
  1 件のコメント
Star Strider
Star Strider 2022 年 11 月 18 日
It would help to have your data. Please post the frequency and gain vectors.

サインインしてコメントする。


Moussa Ihab
Moussa Ihab 2022 年 11 月 18 日
freq=
Columns 1 through 11
2400 2470 2540 2610 2680 2750 2820 2890 2960 3030 3100
Columns 12 through 22
3170 3240 3310 3380 3450 3520 3590 3660 3730 3800 3870
Columns 23 through 33
3940 4010 4080 4150 4220 4290 4360 4430 4500 4570 4640
Columns 34 through 41
4710 4780 4850 4920 4990 5060 5130 5200
gain=
Columns 1 through 13
-18.6680 -14.4060 -10.2174 -6.4374 -3.7861 -2.4066 -1.7975 -1.6368 -1.7193 -1.7336 -1.7873 -1.4274 -1.7189
Columns 14 through 26
-1.3573 -1.6372 -1.3042 -1.2758 -1.1481 -1.0076 -1.1632 -1.0454 -1.1071 -1.1686 -1.1704 -1.3548 -1.2611
Columns 27 through 39
-1.4895 -1.3286 -1.6556 -1.4886 -2.1840 -2.1738 -3.9351 -3.9649 -6.5240 -7.0555 -10.1172 -11.1532 -13.9887
Columns 40 through 41
-16.3625 -18.6674

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by