close all; clear; clc;

%% Real
% Show that for real signals the half-fft is lossless
N = 117;
dt = 0.1; t = (0:(N-1)) * dt; % Time axis
v = linspace(0,1-1/N,N) + (mod(N,2)/N-1)/2; v = v/dt; % Frequency axis
vnq = 1/2/dt; % Nyquist
Np = 5; % Number of plots
Nf = 3; % Number of frequency components

figure(1); clf; cm = lines;
for j = 1:Np
    fs = rand(Nf,1).*vnq; phis = rand(Nf,1).*pi; as = randn(Nf,1);
    ys = as.*sin(2*pi*t.*fs+phis);
    y = sum(ys); % Composite signal
    y = y.*hann(N)'; % Window just to make it look pretty
    Y = fft(y)/N; % Fourier transform
    
    % Half spectrum
    Yh = Y(1:ceil((N+1)/2));
    Yh(2:ceil(N/2)) = Yh(2:ceil(N/2))*2; % Adjust amplitude of the middle samples
    vh = linspace(0, (1-mod(N,2)/N)/2, ceil((N+1)/2)); vh = vh/dt; % Frequency axis
    
    % Reconstruct full spectrum
    Yl = Yh; Yl(2:ceil(N/2)) = Yl(2:ceil(N/2))/2;
    Yr = fliplr(conj(Yl(2:ceil(N/2))));
    Y2 = [Yl,Yr];
    
    % Reconstruct signal
    y2 = real(ifft(Y2*N));
    
    s = subplot(Np,2,2*j-1);
    plot(t,y2,t,y,'.')
    if j~=Np
        s.XTickLabel = '';
    else
        xlabel('time');
    end
    grid on; box on;
    xlim([0,max(t)]);
    
    Sh = abs(Yh); Sh = Sh./max(Sh);
    s = subplot(Np,2,2*j);
    patch([0 vh vnq],[0 Sh 0],cm(1,:),...
        'Edgecolor',cm(1,:),'FaceAlpha',0.2);
    
    s.YTickLabel = '';
    if j~=Np
        s.XTickLabel = '';
    else
        xlabel('frequency');
    end
    grid on; box on;
    ylim([0,1.1]);
    
    drawnow;
end