Find and plot 3dB and 6dB downpoints of spectral peak frequency in a PSD

I have a PSD and find its peak with:

 Pyy, fyy]= periodogram(filterLow_hyd1,[],[],fs);
 [pk, loc] = findpeaks(Pyy,'Npeaks',1,'SortStr','descend');
 Peak_frequency = fyy(loc);
 pkdB = 10 * log10(pk);

How do I find its 3dB/6dB downpoints?

Answer 1

[Star Strider]

This is robust in the event you want to find the -3dB and -6dB points on more than one peak:

fyy = linspace(0, 50, 250);                             % Create Data
Pyy = 0.5*exp(-(fyy-15).^2) + 0.3*exp(-(fyy-40).^2);        % Create Data

[pk,loc] = findpeaks(Pyy,'Npeaks',1,'SortStr','descend');

db3c = 10^(-3/10);                                      % Relative Magnitude At -3dB
db6c = 10^(-6/10);                                      % Relative Magnitude At -6dB
ofst = 10;
for k1 = 1:length(pk)
    varmtx = [Pyy(loc(k1)-ofst:loc(k1)); fyy(loc(k1)-ofst:loc(k1)); Pyy(loc(k1):loc(k1)+ofst); fyy(loc(k1):loc(k1)+ofst)];
    dBpts(k1,1:2) = interp1(varmtx(1,:), varmtx(2,:), pk(k1)*[db6c db3c], 'linear','extrap');
    dBpts(k1,3:4) = interp1(varmtx(3,:), varmtx(4,:), pk(k1)*[db6c db3c], 'linear','extrap');
end

figure(1)
plot(fyy, Pyy)
hold on
for k1 = 1:length(pk)
    plot(dBpts(k1,:), pk(k1)*[db6c db3c db6c db3c], 'r+')
end
hold off
grid

for k1 = 1:length(pk)
    fprintf(1, '\n\t-6dB frequencies = %.3f, %.3f\n', dBpts(k1,[1 3]))
    fprintf(1, '\n\t-3dB frequencies = %.3f, %.3f\n', dBpts(k1,[2 4]))
end

[Butterflyfish]

Dear Star, thank you so much! It is exactly what I would like and your code works with the created data but it doesn't work with my data. I get the error: 'Interpolation requires at least two sample points in each dimension.' I am not sure what it means: it is that my highest peak doesn't have enough points to define the thresholds? I attach here my Pyy and fyy if you are able to have a look. Thanks again for helping!

[Star Strider]

It would have helped to have had your data earlier. My code works with your data without error. You simply have to transpose your data to row vectors:

V1 = load('Butterflyfish Pyy.mat');
V2 = load('Butterflyfish fyy.mat');
Pyy = V1.Pyy';
fyy = V2.fyy';

and change the offset:

ofst = 2;

It produces the plot and this output:

  -6dB frequencies = 59.628, 75.712

  -3dB frequencies = 60.987, 71.882

In order to see the details, you need to add an axis call to the plot, after the rest of the plot code:

axis([0  100    ylim])

[Butterflyfish]

Beauty! Thanks a million, that works perfectly. Sorry for giving the data only afterwards.

[Star Strider]

Thank you! My pleasure!

Your data only necessitated a couple tweaks in my code, one of which (that it needs row vectors) I probably should have mentioned at the outset. The change in ‘ofst’ is the sort of tweak necessary when dealing with real-world data.