%
%   calculation of the spectrum of sinusoid using FFT
%   
%   The input and output time signals, as well as the amplitude
%   spectrum are shown for:
%
%         y = a sin(2pi ft) + b
%
%    call: y = sine_test(a,b,f1)
%
%    default: a=b = 0.2 and f=100 Hz   (dt = 0.0001 s)
%
function y=Chapter3_sine_spectrum(a,b,f1)

if nargin<3 f1=100; end
if nargin<2 b=2; end
if nargin<1 a=0.2; end

dt = 0.0001;

t=[0:dt:2];
fmax = 0.5/dt

x = sin(2*pi*f1*t);
y=a*x+b;

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('f=%3d Hz      [a,b]=[%2.1f, %2.1f]',f1,a,b); 
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);     % the first halve of the spectral representation
fty=fft(y,4096);   fty=fty(1:2048);
ftxa=abs(ftx); ftya=abs(fty);
df = fmax/2048;
freq=[0:df:fmax-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]);
axis([0.001*fmax 0.2*fmax 1 max(ftya)]);

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');

dc=ftya(1);
s=sprintf('dc=%4.1f      ',dc); 
text(400,2, s,'fontsize',11, 'fontweight','bold');

end