Identification of fundamental frequency and not harmonics in a FFT - using peak locations

Hello

I would like to generate my results as similar as I can to the images in red at Figure 6 | Feature extraction of milling chatter based on optimized variational mode decomposition and multi-scale permutation entropy | SpringerLink

I alredy got the FFT from my noisy data and use the command findpeaks I have obtain the following

plot(f1,P1,'--b','LineWidth',1.5) % Frequency domain of the signal

Fz=60; % Fundamental frequency
tooth=Fz;  %Expected peaks separation
hold on

[pkt,lct] =findpeaks(P1,f1,'SortStr','descend','NPeaks',6,'MinPeakDistance',tooth/5)
% pkt is amplitude, and lct is the corresponding frequency
plot(lct,pkt,'r*','MarkerSize',10)

%%
round([pkt,lct'])
ans =
          78         240
          18          60
          14         180
          14         960
          13        1050
          12         480

You can see, only one of these peak is a non harmonics. I want to modify it so:

1) Identify of the group[pkt,lct], which ones are harmonics of the 60Hz component, label each one as harmonics or nor harmonic.

2) Plot harmonics with one symbols and non harmonic with other symbol.

Warning:Some fft will have more than one non harmonics, or can all the peaks be harmonics, so the group data could be something smaller than 6 items.

Attached one image of the current state of the code but as you may notice, it is not working as expected.

Answer 1

[Mathieu NOE]

hello

try this

of course I cannot plot the spectrum as this is not provided here

I just let you see how to do the separation between harmonic and non harmonic data.

BTW , the sorting of the data by findpeaks does not bring any added value here - not needed in fact;

here I start from the data : round([pkt,lct'])

% round([pkt,lct'])
data =[
          78         240
          18          60
          14         180
          14         960
          13        1050
          12         480];
pkt =  data(:,1);
lct =  data(:,2);

% harmonic / non harmonic group segregation
f_base = min(lct);
f_ratio = (lct./f_base);
f_residue = abs(f_ratio-round(f_ratio));

tol = 1e-2; % 1% tolerance
ind_harmo = find(f_residue<=tol);
ind_non_harmo = find(f_residue>tol);

% harmonic group
pkt_harmo =  data(ind_harmo,1);
lct_harmo =  data(ind_harmo,2);

% non harmonic group
pkt_non_harmo =  data(ind_non_harmo,1);
lct_non_harmo =  data(ind_non_harmo,2);

figure(1)
plot(lct_harmo,pkt_harmo,'dr',lct_non_harmo,pkt_non_harmo,'*k');
legend('harmo peaks','non harmo peaks');

[John Navarro]

Thanks so much. It is working great.

The sorting is because later I need to know the location of each peak. In my project it is very relevant to know for example if the peak at 240Hz is the highest or not. Or should I do the sorting in another line, instead that in the findpeaks command?

[Mathieu NOE]

hello

I think you can do the sorting afterwards , but if you need it anyway in your code , where you do it has little importance