fs = 44100;
f0 = 440;

maxharmonic = (fs/2) / f0 - 1;
maxharmonic = 9;

saw_polynomial = chebpolysum([0, 1 ./ (1:maxharmonic)]);

sinusoid = sin(2*pi*f0*(0:(1/fs):3));

decaying_sinusoid = exp(-linspace(5,0,length(sinusoid))) ...
                      .* sinusoid;

xticks = -1:0.01:1;
plot(xticks,polyval(saw_polynomial,xticks));

shapedsound = polyval(saw_polynomial,decaying_sinusoid);
%shapedsound = polyval(saw_polynomial,sinusoid);