function x = mycos(A, f, th)
% MYCOS plots cosine(s) with amplitude(s) A, frequency(s) f,
% and angle(s) th. It also plots Fourier-transform of the signal.
%
%    x = mycos(A, f, th);
%    
%    x  (T x 1)  (sum of) cosine signal(s)
%    A  (N x 1)  amplitude(s)
%    f  (N x 1)  frequency(s)
%    th (N x 1)  angle(s)
%
%    x = SUM_n A_n cos(2 pi f_n t + th_n)
%
% Sampling frequency is set to fs=10000 Hz.
%
% Example
%    % one cosine
%    x  = mycos(3, 200, 0);
%    soundsc(x, 10000);
%    % harmonic overtones
%    y  = mycos([3 2 5 7],[200 400 600 800],[0 pi/6 0.45 -0.2]);
%    soundsc(y, 10000);
%    % non-harmonic
%    z  = mycos(rand(50,1),3000*rand(50,1),2*pi*rand(50,1));
%    soundsc(z, 10000);

% Jukka Parviainen, parvi@hut.fi 2005-11-10

if ((nargin < 3) | (nargin > 3))
  error('Wrong number of arguments. Type "help mycos"');
end;
if (sum(size(A) ~= size(f)) | sum(size(f) ~= size(th)))
  error('Arguments A, f, and th not of same size. Type "help mycos"');
end;
if (prod(size(A)) > 1000)
  error('Too big size of arguments.');
end;

% Initialization 
fs   = 10000;      % sampling frequency
T    = 5000;       % number of samples
len  = T/fs;       % length of signal in seconds
n    = [0 : T-1]'; % index vector
t    = n/fs;       % time vector

% Aliasing?
if (sum(f >= fs/2) > 0)
  warning('One or more components of f over fs/2 - aliasing');
end;

% Signal generation
x    = zeros(T,1);
for i = [1 : length(A(:))]
  tmp = A(i) * cos(2*pi*f(i)*t + th(i));
  x   = x + tmp;
end;

% Plotting signal
%MATLAB: clf;              % clear figure
subplot(211);     % upper graph
plot(t, x);
grid on;
xlabel('Time (sec)');
ylabel('Amplitude');
hold on;
plot(t, x, 'ro');
%MATLAB: set(gca,'XLim',[0 0.05]);     % zoom only to 0 .. 0.05 s
title('Signal in time domain');

% Fourier-transform
xF   = fft(x);          % Discrete Fourier Transform (DFT)
w    = fs*[0 : T-1]/T;  % frequency axis

% Spectrum
subplot(212);
plot(w, abs(xF));       % linear scale
hold on;
plot(w, abs(xF), 'ro');
%plot(w, 20*log10(abs(xF)));       % logaritmic scale
%hold on;
%plot(w, 20*log10(abs(xF)), 'ro');
axis([0 fs/2 0 1.1*max(abs(xF))]);
%set(gca, 'XLim', [0 fs/2]);       % logaritmic scale
grid on;
xlabel('Frequency (Hz)');
ylabel('Amplitude');
title('Spectrum in frequency domain');
