Extract maxima of a fitted curve

%%% Parameters %%%

PPmax = 100;   %Max # of points per period
PPmin = 2;   %Min # of points per period
incs = 2;   %number of samples increment
F1 = 1;   %frequency 1
F2 = 1.2;   %frequency 2
NP = 04;   %number of periods of shorter frequency
Nz = 1000;   % number of zeros to pad with
o = 0;   %phase shift

                                
N = PPmax*NP;   % max number of sample points
T = abs(1/(F1-F2));   %one period of shorter frequency
L = NP*T;   %length of signal   
dt = L/(N-1);   %time step
f_basis = dt*(Nz+N-1);   %used for converting to frequency
f_axis = (0:1:N+Nz-1)/(f_basis);   % frequency axis for FFT    
t = 0:dt:L; 

%%% Signal %%%

y = cos(t.*pi.*(F1+F2)+o).*cos(t.*pi.*(F1-F2)+o);   

%%% Windowing %%%

win = gausswin(N);   
y = y'.*win;

%%% Zero padding %%%

y = [y' zeros(1,Nz)];

%%% Fourier Transform %%%

Y = (abs(fft(y))).^2;

%%% Fitting %%%
    
[xData, yData] = prepareCurveData( U_axis, U );

% Set up fittype and options.
ft = fittype( 'smoothingspline' );
opts = fitoptions( 'Method', 'SmoothingSpline' );
opts.SmoothingParam = 1;

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
    
%%% Extract Maxima %%%

yfitted = feval(fitresult,xData);   
[ypk,idx] = findpeaks(yfitted);
[~,loci] = maxk(ypk,2);
loc = idx(loci); 
F1_found = (min(loc))/f_basis_loop
F2_found = (max(loc))/f_basis_loop

plot(fitresult)

Answer 1

[Mathieu NOE]

hello

your code cannot be run because you did not provide prepareCurveData function and the data U_axis, U

[xData, yData] = prepareCurveData( U_axis, U );

also I don't understand why you need to fit these data

finding and plotting the two frequencies peaks is indeed quite simple here , and your code could be made a little more concise.

%%% Parameters %%%
PPmax = 100;   %Max # of points per period
PPmin = 2;   %Min # of points per period
incs = 2;   %number of samples increment
F1 = 1;   %frequency 1
F2 = 1.2;   %frequency 2
NP = 04;   %number of periods of shorter frequency
Nz = 1000;   % number of zeros to pad with
o = 0;   %phase shift
                                
N = PPmax*NP;   % max number of sample points
T = abs(1/(F1-F2));   %one period of shorter frequency
L = NP*T;   %length of signal   
dt = L/(N-1);   %time step
f_basis = dt*(Nz+N-1);   %used for converting to frequency
f_axis = (0:1:N+Nz-1)/(f_basis);   % frequency axis for FFT    
t = 0:dt:L; 
%%% Signal %%%
y = cos(t.*pi.*(F1+F2)+o).*cos(t.*pi.*(F1-F2)+o);   
%%% Windowing %%%
win = gausswin(N);   
y = y'.*win;
%%% Zero padding %%%
y = [y' zeros(1,Nz)];
%%% Fourier Transform %%%
Y = (abs(fft(y))).^2;
m = length(f_axis);

%% consider only first half freq range (second half symmetrical)
f_axis = f_axis(1:m/2);
Y = Y(1:m/2);

[ypk,idx] = findpeaks(Y,'MinPeakHeight',max(Y)/4);
F1_found = f_axis(min(idx))
F2_found = f_axis(max(idx))
plot(f_axis,Y,f_axis(idx),Y(idx),'dr')

[Jacob Duryee-Feeney]

I need the spline because sometimes the true peak is between two of the data points of the FFT. It helps increase the accuracy a bit.

[Mathieu NOE]

ok, I implemented your idea , but with the "regular" spline function (as I don't have the curve fit tbx)

yes we can see that the second computation is slightly closer to the true value

F1_found = 0.998213009292352

F2_found = 1.197855611150822

F1_found = 1.000511197401833

F2_found = 1.201852460036877

clc;
close all;
clear all;
%%% Parameters %%%
PPmax = 100;   %Max # of points per period
PPmin = 2;   %Min # of points per period
incs = 2;   %number of samples increment
F1 = 1;   %frequency 1
F2 = 1.2;   %frequency 2
NP = 04;   %number of periods of shorter frequency
Nz = 1000;   % number of zeros to pad with
o = 0;   %phase shift
                                
N = PPmax*NP;   % max number of sample points
T = abs(1/(F1-F2));   %one period of shorter frequency
L = NP*T;   %length of signal   
dt = L/(N-1);   %time step
f_basis = dt*(Nz+N-1);   %used for converting to frequency
f_axis = (0:1:N+Nz-1)/(f_basis);   % frequency axis for FFT    
t = 0:dt:L; 
%%% Signal %%%
y = cos(t.*pi.*(F1+F2)+o).*cos(t.*pi.*(F1-F2)+o);   
%%% Windowing %%%
win = gausswin(N);   
y = y'.*win;
%%% Zero padding %%%
y = [y' zeros(1,Nz)];
%%% Fourier Transform %%%
Y = (abs(fft(y))).^2;
m = length(f_axis);

%% consider only first half freq range (second half symmetrical)
f_axis = f_axis(1:m/2);
Y = Y(1:m/2);
[ypk,idx] = findpeaks(Y,'MinPeakHeight',max(Y)/4);
F1_found = f_axis(min(idx))
F2_found = f_axis(max(idx))

% refine frequency range around F1 and F2 first computation
f_min = 0.5*F1_found;
f_max = 1.5*F2_found;

ind = find(f_axis>=f_min & f_axis<=f_max);
f_axis = f_axis(ind);
Y = Y(ind);
f_axis_refine = linspace(f_min,f_max,1000);
% spline smooth 
Yq = spline(f_axis,Y,f_axis_refine);

[ypk,idx] = findpeaks(Yq,'MinPeakHeight',max(Yq)/4);
F1_found = f_axis_refine(min(idx))
F2_found = f_axis_refine(max(idx))
plot(f_axis,Y,'-*k',f_axis_refine,Yq,'-g',f_axis_refine(idx),Yq(idx),'dr')
legend('raw ftt','spline smooth','peaks');