function c = sinesamples(fs, f1, f2, time, pers, apers)
%c = sinesamples(fs, f1, f2, time, pers, apers)
%makes a variety of sine based samples at sampling frequency fs.
%for samples that involve 2 frequencies will go between f1 and f2.  f1 < f2
%each sample will be 'time' seconds long
%pers is the number of periods the fm goes through
%apers if the number of periods the am goes through
%returns samples in cell array c

%this script will synthesize a couple of fm sweeps (including static tones) and 
%process them with varying parameters to see how they sound.

x1 = createsine(880, fs, time);

%create a linear fm sweep from high A (880) to B
s = 1.059463; %half step ratio
A = f1;
B = f2;

len = time;
t = (0:(1/fs):len)';
x2 = sin(2*pi*(A*t+(B-A)/(len*2)*t.^2));

%create alternating tones between A and B

len = time;
t = (0:(1/fs):(len))';
lenT = length(t);
x3 = zeros(lenT,1);
p = 0;

fmTsamps = length(t)/pers;

for i = 1:length(t)
    x3(i) = sin(p);
    if mod(i, fmTsamps) <= fmTsamps/2
        p = p + 2*pi/(fs/A);
    else
        p = p + 2*pi/(fs/B);
    end
end

%create linear sweeps between A and B

%first attempt, does not preserve phase when changing sweeping rate
% lam = lambda * ones(length(t)/4,1);
% i = 1:(length(t)/4);
% x31 = sin(2.*pi.*(A.*t(i)+(lam/2).*t(i).^2));
% i = (length(t)/4 + 1):(length(t)/2);
% x32 = sin(2.*pi.*(A.*t(i)-(lam/2).*t(i).^2 + lam*len/2.*t(i)));
% i = (length(t)/2 + 1):(3*length(t)/4);
% x33 = sin(2.*pi.*(A.*t(i)+(lam/2).*t(i).^2 - lam*len/2.*t(i)));
% i = (3*length(t)/4 + 1):length(t);
% x34 = sin(2.*pi.*(A.*t(i)-(lam/2).*t(i).^2 + lam*len.*t(i)));
% x3 = [x31;x32;x33;x34];

len = time;
t = (0:(1/fs):(len))';
lenT = length(t);
x4 = zeros(lenT,1);
p = 0;
%define dp, this is how much we change the phase by each time;
dp = 0;
ddp = 2*pi*4*(B-A)/(len*fs^2);
for i = 1:length(t)
    x4(i) = sin(p);
    p = p + 2*pi*A/fs + dp;
    if mod(i, fmTsamps) <= fmTsamps/2
        dp = dp + ddp;
    else
        dp = dp - ddp;
    end
end

%create a sine fm sweep between high A (880) and B

len = time;
t = (0:(1/fs):len)';

mid = (A + B)/2;
amp = (B - A)/2;
w = pers*2*pi/len;
x5 = sin(2.*pi.*mid.*t + amp*sin(w.*t-pi/2));

c = cell(5, 6);
c{1, 1} = x1;
c{2, 1} = x2;
c{3, 1} = x3;
c{4, 1} = x4;
c{5, 1} = x5;

%now any of these can be amplitude modulated
%lets do 5 modulations:
%linear up
%linear down
%exp up
%exp down
%sine

%but the expontials grow too linear, so lets use a sigmoid

%linear up
%still have t
env = linspace(0, 1, length(t))';
c{1, 2} = x1 .* env;
c{2, 2} = x2 .* env;
c{3, 2} = x3 .* env;
c{4, 2} = x4 .* env;
c{5, 2} = x5 .* env;

%linear down
env = linspace(1,0,length(t))';
c{1, 3} = x1 .* env;
c{2, 3} = x2 .* env;
c{3, 3} = x3 .* env;
c{4, 3} = x4 .* env;
c{5, 3} = x5 .* env;

%sigmoid up
y = linspace(-6,6,length(t))';
env = 1./(1+exp(-y));
env = env / max(env);
c{1, 4} = x1 .* env;
c{2, 4} = x2 .* env;
c{3, 4} = x3 .* env;
c{4, 4} = x4 .* env;
c{5, 4} = x5 .* env;

%exponential down
env = flipud(env);
c{1, 5} = x1 .* env;
c{2, 5} = x2 .* env;
c{3, 5} = x3 .* env;
c{4, 5} = x4 .* env;
c{5, 5} = x5 .* env;

%sine modulation.  Lets do 2 periods
y = linspace(1,3,length(t))';
env = abs(sin(apers*2*pi/time*y));
c{1, 6} = x1 .* env;
c{2, 6} = x2 .* env;
c{3, 6} = x3 .* env;
c{4, 6} = x4 .* env;
c{5, 6} = x5 .* env;