% plot spectrum for linear and nonlinear time series
% example for nonlinear dynamics and chaos course.  Eli Tziperman
clear all;

N=1000;
a=3;
b=10;
omega1=91*(2*pi/N);
omega2=25*(2*pi/N);

for i=1:N
  % change the sin() to sin()^3 or sin()^5 to see the effects of
  % nonlinearity on the spectrum;  change the ratio of a and b to
  % see the relation between the amplitude of a given frequency in the
  % time series and its amplitude in the spectrum:
  x(i)=a*sin(omega1*i)+b*sin(omega2*i);
end

% plot the time series:
figure(1)
plot(x)
xlabel('time')
ylabel('x(t)')
title('time series')

figure(2)
% Calculate the spectrum if the spectrum command is not available:
% First, calculate the Fourier transform:
Y=fft(x);
Y(1)=[];
n=2*floor(length(Y)/2);
% Second, calculate the square of the absolute value of the Foureir
% coefficients:
power=abs(Y(1:n/2)).^2;
nyquist=1/2;
% Third, calculate the frequencies:
freq=(1:n/2)/(n/2)*nyquist;

% calc the power spectrum using the spectrum command, if the right
% matlab toolbox is available: 
% P=spectrum(x); Power=P(:,1)

% plot the spectrum:
% multiply freq by 2*pi to put in rad/s
semilogy(freq*2*pi,power)
xlabel('frequency (rad/s)')
ylabel('Power Spectrum')
%title('power spectrum (roughly squared fourier coefficients) on a log scale')
title(['a=' num2str(a) ', \omega_1=' num2str(omega1)...
       ', b=' num2str(b) ', \omega_2=' num2str(omega2)])
