f = 2000;
fs = 44100;
fs_upsampled = 176400;
s = 1:(fs * 2);
s_upsampled = 1:(fs_upsampled * 2);

% naive
naive = 1 - 2 * mod(s*f/fs, 1);

% ideal
ideal = naive * 0;
f1 = f;
i = 0;
while f1 < fs / 2
    i = i + 1;
    f1 = f * i;
    ideal = ideal + sin(2 * pi * s * f1 / fs) / i;
end
ideal = ideal / norm(ideal) * norm(naive);

% naive + oversampling
naive_upsampled = 1 - 2 * mod(s_upsampled*f/fs_upsampled, 1);
naive_upsampled = resample(naive_upsampled, 1, 4);

figure;
hold on;
plot(ideal, 'r')
plot(naive_upsampled, 'g')
plot(naive)

figure;
hold on;
plot(abs(fft(naive)))
plot(abs(fft(naive_upsampled)), 'g')
plot(abs(fft(ideal)), 'r')


sound(ideal, fs)
sound(naive, fs)
sound(naive_upsampled, fs)