Extract maxima of a fitted curve

3 views (last 30 days)
%%% 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)

Accepted Answer

Mathieu NOE
Mathieu NOE on 26 Apr 2022
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')
  2 Comments
Mathieu NOE
Mathieu NOE on 26 Apr 2022
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');

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!