close all; clear; clc;

%%
y = @(x) cos(2*pi*2*x); % 2Hz signal

Ns = unique(round(10.^linspace(1.5,3.5,200)));
t_fast = NaN(1,numel(Ns)); d_fast = NaN(1,numel(Ns));
t_naive = NaN(1,numel(Ns)); d_naive = NaN(1,numel(Ns));
%t_fast = zeros(1,numel(Ns)); t_naive = zeros(1,numel(Ns)); d_fast = zeros(1,numel(Ns)); d_naive = zeros(1,numel(Ns));
for Ni = 1:numel(Ns)
    N = Ns(Ni);
    t = linspace(0,5,N+1); t(end) = []; % Time axis
    
    % Fast
    tic
    for j=1:10
        fast = fft(y(t));
    end
    t_fast(Ni) = toc*1e3/10;
    d_fast(Ni) = sqrt(median(imag(fast).^2)); % Root median square imag
    
    % Naive
    dt = t(2)-t(1);
    v = linspace(0,1-1/N,N) + (mod(N,2)/N-1)/2; v = v/dt;
    tic
    naive = sum(bsxfun(@times,y(t),exp(2 * pi * 1i .* v' .* (t-t(1)))),2)';
    t_naive(Ni) = toc*1e3;
    d_naive(Ni) = sqrt(median(imag(naive).^2));
    if ~mod(Ni,round(numel(Ns)/20)), fprintf('%.1f%% done\n',Ni/numel(Ns)*100); end
end

%% Plot
% Fit with lines on log-log scale
p1 = polyfit(log10(Ns),log10(d_naive),1); p2 = polyfit(log10(Ns),log10(d_fast),1);

figure(1); clf;
subplot(1,2,1);
h = loglog(Ns,d_naive,Ns,d_fast); hold on;
%loglog(Ns,exp(p1(2)).*Ns.^p1(1),'k--',Ns,exp(p2(2)).*Ns.^p2(1),'k--');
loglog(Ns,10.^(p1(2)).*Ns.^p1(1),'k--',Ns,10.^(p2(2)).*Ns.^p2(1),'k--');
box on; xlim([min(Ns) max(Ns)]); axis square;
xlabel('N'); ylabel('rms imag(FT(x))');
%title(sprintf('naive ~ O(N^{%.1f}), fast ~ O(N^{%.1f})',p1(1),p2(1)));
legend(h,'naive','fast','location','northwest');
grid on;

%
s1 = t_naive(Ns>100) / ( Ns(Ns>100).^2 );
s2 = t_fast(Ns>100) / ( Ns(Ns>100) .* log2(Ns(Ns>100)) );

subplot(1,2,2);
h = loglog(Ns,t_naive,Ns,t_fast); hold on;
loglog(Ns,s1*Ns.^2,'k--',Ns,s2*Ns.*log2(Ns),'k--');
box on; xlim([min(Ns) max(Ns)]); axis square;
xlabel('N'); ylabel('execution time / ms');
%legend(h,'naive','fast','location','northwest');
grid on;
