clear
close all
load clipping_measurement.mat

T= N/fs;    %sweep duration 
f1=20;  
f2=20e3;
Rate=log(f2/f1);    %sweep rate
IR = fftshift(ifft(fft(Vout,2*N-1)./fft(Vin,2*N-1)));
%%
N = 10;
dt = T*log( (1:N) )./Rate *fs;
pks = round(length(IR)/2-dt);
%%
Csin = zeros(N);
for n = 0 : N        % Sine
    if rem(n,2) == 0   % Even
        k = (0:n/2-1);
        c0 = factorial(n)./(factorial(n-k).*factorial(k));        
        k0 = (n-2*k);
        Csin(n+1, k0+1) =  2/2^n*c0.*(-1).^(n/2-k)*1i; 
        Csin(n+1,1) = Csin(n+1,1) + 1/(2.^n)*nchoosek(n,n/2) ;
    elseif rem(n,2) == 1 % Odd
        k = (0:(n-1)/2);
        c0 = factorial(n)./(factorial(n-k).*factorial(k));        
        k0 = (n-2*k);
        Csin(n+1, k0+1) = 2/2^n*c0.*(-1).^((n-1)/2-k);
    end
end
%%
A = Csin(2:end,2:end).';
L = 256*8;
h = zeros(L,N);
for i=1:N
   h(:,i)= flipud(IR(pks(i)-L/2 : pks(i)+L/2-1 ));
end
%%
H = fft(h,[],1);
H_volterra = (pinv(A)*H.').';
h_volterra = ifft(H_volterra,L,1,'symmetric');

t = (0:fs-1)'/fs;
f0 = 500;
in = sin(2*pi*f0*t);

out_volterra = zeros(length(in),1);
for n = 1 : N
    out_volterra = out_volterra + conv( in.^n, h_volterra(:,n),'same' );
end
out_volterra = out_volterra(1:length(in));

figure
subplot(1,2,1)
plot(t,out_volterra,'LineWidth',1.5)
grid on
xlim([10/f0,12/f0])
xlabel('t -> [s]')

subplot(1,2,2)
f = (0:length(in)-1)/length(in)*fs;
hold on 
plot(f,(abs(fft(out_volterra))),'LineWidth',1.5)
grid on
xlim([0,f0*10])
xlabel('f -> [Hz]')
