%
%   simulation of the spectrum of a third-order nonlinear system
%   in response to two sine waves of different frequencies.
%   The input and output time signals, as well as the amplitude
%   spectra are determined.
%
%         y = ax + bx^2 + cx^3
%
%   call: y = nonlinearsines(a,b,c,f1,f2)
%    default: a=b=c=0.2 and f1=100 Hz and f2=175 Hz;
%
function y=nonlinearsines(a,b,c,f1,f2)
if nargin<5 f2=175; end
if nargin<4 f1=100; end
if nargin<3 c=0.2; end
if nargin<2 b=0.2; end
if nargin<2 a=0.2; end

t=[0:0.0001:2];
x1 = sin(2*pi*f1*t);
x2 = sin(2*pi*f2*t);
x=x1+x2;

y=a*x+b*x.^2+c*x.^3;

figure(1)
clf
plot(t,x,'r-');
hold on
plot(t,y,'k-');
axis([0 0.1 -3 3]);
set(gca,'fontsize',12,'fontweight','bold');
xlabel('Time (s)', 'fontsize',13,'fontweight','bold');
ylabel('Amplitude (a.u.)', 'fontsize',13,'fontweight','bold');

s=sprintf('f1=%3d Hz, f2=%3d Hz      [a,b,c]=[%2.1f, %2.1f,%2.1f]',f1,f2,a,b,c); 
text(0.01,-2.5, s,'fontsize',11, 'fontweight','bold');
print -djpeg fig4-1a

%
% determine the input and output spectra ...
%
ftx=fft(x,4096);   ftx=ftx(1:2048);
fty=fft(y,4096);   fty=fty(1:2048);
ftxa=abs(ftx); ftya=abs(fty);
df = 5000/2048;
freq=[0:df:5000-df];

figure(2); clf;
%semilogx(freq, ftxa,'r-');
loglog(freq, ftxa,'r-');
hold on
%semilogx(freq, ftya,'k-');
loglog(freq, ftya,'k-');
axis([10 1000 1 2000]);

set(gca,'fontsize',12,'fontweight','bold');
xlabel('Frequency (Hz)', 'fontsize',13,'fontweight','bold');
ylabel('Log-Amplitude', 'fontsize',13,'fontweight','bold');
plot([15,25],[5,5],'r-','linewidth',2); text(27,5,'  Input spectrum','fontsize',11,'fontweight','bold');
plot([15,25],[4,4],'k-','linewidth',2); text(27,3.5,'Output spectrum','fontsize',11,'fontweight','bold');

print -djpeg fig4-1b


end