function s = afwosc(P,F,A,mr,or)
%
% Amplitude- and Frequency-modulated signals based on Wavetable Oscillator
%
% s = wosc(P,F,A,mr,or)
%
%  P = one period of a sampled signal
%  F = vector of output frequencies
%  A = vector of output amplitudes
%  mr = sampling rate of modulation signals (e.g. 200)
%  or = sampling rate of output signal (e.g. 11025)
%

lP = length(P);
lF = length(F); lA = length(A);  % should check that these are the same
P1 = [P P(1)];  % in case of interpolation across the wrap
%
n  = floor(or*(lF/mr));  % number of output samples
%
% Index into modulation vectors per output sample:
iM = floor(mr*(0:(n-1))/or)+1;
% Modulation values per output sample:
mF = F(iM);  
mA = A(iM);
%
% Added: up- and down-ramps to amplitude modulation,
%   to avoid clicks at edges of pitched regions
Ax = [A(2:length(A)); -1];
uA = find((A==0)&(Ax > 0)) - 1 ;
dA = find((A > 0)& (Ax == 0)) - 1;
Ramps = ones(length(mF),1);
ups = round(.02*or); upramp = (1:ups)/ups;
downs = round(.03*or); downramp = (downs:-1:1)/downs;
spf = or/mr;

for n=1:length(uA)
    i1 = floor(uA(n)*spf);
    Ramps(i1:(i1+ups-1)) = upramp;
end

for n=1:length(dA)
    i1 = ceil((dA(n)+1)*spf)-downs+1;
    Ramps(i1:(i1+downs-1)) = downramp;
end

% Wavetable increments per output sample:

inc = lP*mF/or;  % increment in P per output sample:
                 % each output sample is 1/or long;
                 % each sample in P is nominally (1/F)/lP long;
                 % thus if intersample distance in P is 1,
                 % intersample distance in output is (1/or)/((1/F)/lP)
% Cumulative increments, modulo wavetable length
%   (with 1 added to deal with matlab index issue)
Cinc = mod(cumsum(inc), lP) + 1;
CincR = round(Cinc);
% Interpolation -- not really necessay, but if you want:
%Cinc1 = floor(Cinc);
%Cinc2 = ceil(Cinc);
%frac2 = Cinc-Cinc1;
%frac1 = 1.0 - frac2;
% s = (frac1*P1(Cinc1)+frac2*P1(Cinc2)).*mA'.*Ramps';
s = P1(CincR).*mA'.*Ramps';
return

