% Nonlinear circuit
%
% -----------------------------------------------------------------------------------
%
% (C) 2017, Krzysztof Sozanski, Digital Signal Processing in Power Electronics Control 
% Circuits, second edition, Springer-Verlag, London, 2017. 
% 
% -----------------------------------------------------------------------------------
%
clear all; close all;
  roz_fon=18; grub_lin=2;
  N=2^12; % number of samples
  fs=12800; % sampling frequency
  f1=50; f_1k=round(f1/(fs/N))*fs/N; % line frequency
  A1=1; 
  t=(0:N-1)/fs; % time vector
  a=[1 0 -0.2 0 0.15 0 -0.11 0 -0.05]; % coefficients
  %% ----- Response for DC -----
  x=((0:N-1)/N)*2-1; % input signal -1...1
  y=a(1)*x+a(2)*x.^2+a(3)*x.^3+a(4)*x.^4+a(5)*x.^5 + ...
    +a(6)*x.^6+a(7)*x.^7+a(8)*x.^8+a(9)*x.^9;
  plot(x,x,x,y,'LineWidth',grub_lin); 
  	set(gca,'FontSize',[roz_fon],'FontWeight','d'),
   	xlabel('x'); ylabel('y, x'); grid on;
  %% ----- Response for sinusoidal signal 100 Hz -----
  x=A1*sin(2*pi*f_1k*t); % input sinusoidal signal 
  y=a(1)*x+a(2)*x.^2+a(3)*x.^3+a(4)*x.^4+a(5)*x.^5+... 
    +a(6)*x.^6+a(7)*x.^7+a(8)*x.^8+a(9)*x.^9;
  % spectrum
  w_y=fft(y)/N*2; w_y_dB=20*log10(abs(w_y)+eps);
  f=(0:N-1)*fs/N; 
  subplot(211);
	plot(t,x,t,y,'LineWidth',grub_lin); grid on;
    set(gca,'FontSize',[roz_fon],'FontWeight','n'),
    set(gca,'Xlim',[0 2/f_1k]); xlabel('Time [s]');
    ylabel('Amplitude'); title('(a)');
  subplot(212);
 	plot(f,w_y_dB,'r','LineWidth',grub_lin+1); grid on;
 	set(gca,'Ylim',[-80 0]);set(gca,'Xlim',[0 11*f_1k]);
 	set(gca,'FontSize',[roz_fon],'FontWeight','n'), 
 	xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]'); 
 	title('(b)'); 
