%
% ssyn.m: 正弦波の生成・合成 (09/11/13; 12/01/26, 19/11/08 YH)
%
%   [y Fs Ndata] = ssyn(f, tlen [, options])
%       f   : 周波数
%       tlen: 時間長（秒）
%       他のパラメタはキーワードパラメタとして指定する。
%       f がベクトルのときは、それらの周波数の正弦波を合成する。
%       キーワードパラメタ
%           'Fs': サンプリングレート（デフォールト 44100）
%           'A' : 振幅（ベクトルの場合、周波数ベクトルの各成分の振幅）
%           'ph': 位相（ベクトルの場合、周波数ベクトルの各成分の位相）
%           'play': （パラメタなし）：音を鳴らす。

function [y, Fs, Ndata] = ssyn(f, tlen, varargin)
%
% 以下、デフォールト値の設定
%
Fs = 44100;                 % サンプリングレート
A = 1.0;                    % 振幅
ph = 0.0;                   % 位相
play = 0;                   % 演奏（デフォールトは鳴らさない）
if (nargin < 2)             % デフォールト時間長：１秒
    tlen = 1.0;
end

i = 1;
vl = length(varargin);
while (i <= vl)
    switch (varargin{i})    % 注意：エラーチェックはしていない
        case {'Fs', 'fs', 'FS'}                             % サンプリングレート
            i = i + 1; Fs = varargin{i};    
        case {'A', 'Amp', 'amp', 'Amplitude', 'amplitude'}  % 振幅
            i = i + 1; A = varargin{i};
        case {'ph', 'Ph', 'Phase', 'phase'}                 % 位相
            i = i + 1; ph = varargin{i};
        case {'Play', 'play'}                               % 音を鳴らす
            play = 1-play;
        otherwise
            error('error in argument %d\n', i+2);
    end
    i = i + 1;
end

dt = 1/Fs;                                  % １サンプルの時間長
if (mod(tlen, dt) ~= 0)
    t = 0:dt:tlen;                          % サンプル時間列
else
    t = 0:dt:(tlen-dt);                     % tlen が dt で割り切れる場合には最後の１点を除く
end
Ndata = length(t);

y = A(1) * sin(2*pi*f(1).*t + ph(1));       % 最初の成分
if (length(f) > 1)
    if (length(ph) < length(f))             % 周波数ベクトルに対応した位相ベクトルの設定
        for i = length(ph)+1:length(f)
            ph(i) = ph(1);
        end
    end
    if (length(A) < length(f))             % 周波数ベクトルに対応した振幅ベクトルの設定
        for i = length(A)+1:length(f)
            A(i) = A(1);
        end
    end
    for i = 2:length(f)                    % ２番目以降の成分の合成
        y = y + A(i) * sin(2*pi*f(i).*t + ph(i));
    end
end

% y = normalize(y);                         % 正規化：以下でインライン処理
ymax = max(abs(y));
if (ymax > 1); y = y / ymax; end

if (play ~= 0); sound(y, Fs); end        % 音の演奏