function w=fresolve(x,K,w_GT)
% USAGE : w=fresolve(x,K,[w_GT]);
% FUNCTION : Super-resolves the largest K peaks of the DFT of the signal x
% using the formula published in a "Tips & Tricks" column of the IEEE
% Signal Processing Magazine (2023). If the ground-truth frequencies, w_GT,
% are given as input, comparisons with the super-resolved frequencies are
% shown as well, together with some visualization.
% If the DFT has less than K peaks, a warning occurs.
%
% INPUT : * x, the samples of the signal
% * K, the number of peaks to resolve
% * w_GT, the ground-truth frequencies (optional), only needed
% for comparison with the super-resolved frequencies and plots
% OUPUT : * w, the (cyclic) frequency after super-resolution (in ]-pi,pi])
%
% DATE : October 2023
% AUTHOR : Thierry Blu, mailto:thierry.blu@m4x.org
% REFERENCE: Guo, R. & Blu, T.,"Super-Resolving a Frequency Band",
% IEEE Signal Processing Magazine, 2023. To appear.
%__________________________________________________________________________
% EXAMPLE OF USE
% % Generation of N samples of a signal made of K frequencies
% K=5;N=30;
% % First, choosing frequencies (reasonably separated from each other)
% sep=5;
% w_GT=sort(2*pi*rand(sep*K,1));
% w_GT=pi-w_GT(1:sep:end);
% % Second, choosing amplitudes (with reasonably close absolute values)
% a_GT=randn(K,2)*[1;1i];
% a_GT=a_GT./abs(a_GT).*(1+0.1*rand(K,1));
% % Third, building the samples of the signal
% x=a_GT.'*exp(1i*w_GT*(0:(N-1)));
%
% % Finally, super-resolving the K DFT peaks
% w=fresolve(x,K,w_GT);
x=x(:);
N=length(x);
n=0:(N-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% Main calculation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%__________________________________________________________________________
% Functions
arg=@(theta)pi-mod(pi-theta,2*pi); % Wrap frequencies in ]-pi,pi]
Fourier=@(w)exp(-i*w*n)*x; % DTFT coefficient at w
trick=@(w)w+2*atan(tan(pi/(2*N))*... % Super-resolution formula
(abs(Fourier(w+pi/N))-abs(Fourier(w-pi/N)))./...
(abs(Fourier(w+pi/N))+abs(Fourier(w-pi/N))));
%__________________________________________________________________________
% Peaks of the DFT of the signal x
y=abs(fft(x));
y=[y(N);y;y(1)]; % Periodic extension of the samples
[ym,j]=findpeaks(y,'SortStr','descend','NPeaks',K); % Peaks found exclude first and last sample
w_DFT=arg(2*pi*mod(j-2,N)/N); % Frequencies of the DFT peaks
% Has findpeaks found the requested K peaks, or less?
if length(j)K_SR))~=0); % from w and from w_GT are consecutive
err=abs(diff([w_all(end)-2*pi;w_all])); % Error between consecutive frequencies
[~,k]=min(err(idx_trans)); % Index of the minimal error when the consecutive
k=idx_trans(k); % frequencies are from w and w_GT
if j(k)>K_SR
j_GT=[j_GT;j(k)-K_SR];
j_SR=[j_SR;j(mod(k-2,length(j))+1)];
else
j_GT=[j_GT;j(mod(k-2,length(j))+1)-K_SR];
j_SR=[j_SR;j(k)];
end
w_all=w_all(setdiff(1:end,... %
[mod(k-2,length(j))+1,k])); % Remove paired frequencies from w_all
j=j(setdiff(1:end,... %
[mod(k-2,length(j))+1,k])); % and their indices in w and w_GT
stop=(max(j)<=K_SR|min(j)>K_SR); % Only remaining frequencies are either SR, or GT
end
if max(j)>K_SR
w_SR=[w(j_SR);NaN(K_GT-length(j_SR),1)];
w_GT=[w_GT(j_GT);w_GT(j-K_SR)];
else
w_SR=[w(j_SR);w(j)];
w_GT=[w_GT(j_GT);NaN(K_SR-length(j_GT),1)];
end
err=arg(w_SR-w_GT)/(2*pi/N); % 2*pi/N normalized error between the SR and the GT frequencies
disp(' ')
matlab_version=version('-release');
if str2num(matlab_version(1:4))<=2019 % variable names are more restricted in older Matlab versions
disp(table(w_SR,w_GT,err,...
'variablenames',{'w_SR' 'w_GT' 'error_x_N/2\pi'}))
else
disp(table(w_SR,w_GT,err,...
'variablenames',{'w_SR' 'w_GT' 'error x N/2\pi'}))
end
disp(' ')
disp('w_SR = super-resolved frequency')
disp('w_GT = ground-truth frequency')
disp(['DFT bin size = ' num2str(2*pi/N) ' (2\pi/N)'])
end