Aliasing Square Wave
My guess is there's a problem with your wrapPhase function.

If I use the R implementation of it below (C8 shown), then I get all sorts of funny results.  That's the second image.  Removing it from the iteration (but including it in the calculation), makes for much cleaner signals (the first image).







R Code Below

#30582

wrapPhase &lt;- function(phi)
{
  return(phi %% 2*pi)
}

# https://en.wikipedia.org/wiki/Piano_key_frequencies
cs &lt;- c(32.7032, 65.4064, 130.813, 261.626, 523.251, 1046.50, 2093.00, 4186.01)
c1 &lt;- cs[1]
c8 &lt;- cs[8]

fs &lt;- 44100

phs &lt;- rep(0,T + 1)
phs_incr &lt;- c8 / fs * 2 * pi

T &lt;- 1000

sample &lt;- rep(0,T)
sqr_wave &lt;- rep(0,T)

for (t in seq(1,T))
{  
  sample[t] &lt;- sin(wrapPhase(phs[t]))

  if (sample[t] &gt; 0)
  {
    sqr_wave[t] &lt;- 1
  }
  else 
  {
      sqr_wave[t] = -1;
  }

  print(sample[t])

  phs[t+1] &lt;- phs[t] +  phs_incr
  #phs[t+1] &lt;- wrapPhase(phs[t+1]);
}

phi &lt;- seq(-pi, pi, 0.001*pi)
phi_wrap &lt;- rep(0, length(phi))

for (k in 1:length(phi))
{
  phi_wrap[k] &lt;- wrapPhase(phi[k])
}

plot(sample[1:50], col="blue", type="l")
lines(sqr_wave[1:50], col="green")