CMA Equaliser and FSK
The confusion comes from the fact that what is tagged as "transmitted" isn't the real transmitted waveform but its baseband representation, which are $\left\{+1,-1\right\}$ symbols in 1 dimension. For 2-FSK, that representation would be $\left\{\left(+1,0\right), \left(0,+1\right)\right\}$ symbols in 2 dimensions.

A common convenient way to handle such 2 dimensional symbols is to represent them with complex number, mapping the coordinate $\left(+1,0\right)$ to the purely real value $+1$ and the coordinate $\left(0,+1\right)$ to the purely imaginary value $j=\sqrt{-1}$. Using that mapping, the symbol generation code for 2-FSK becomes:

if strcmp(tx_type, "fsk")
  tx_bits = round(rand(1,N));
  s0 = zeros(N,1);
  s0M = zeros(N*M,1);

  for i=1:N
    if tx_bits(i)
      s0(i) = 1i;
    else
      s0(i) = 1;
    end
  end

  k = 1;
  for i=1:M:N*M
   s0M(i:i+M-1) = s0(k);
    k ++;
  end
end


The problem then is that the CMA weight update rule from the original post is given for a purely real signal. The adjusted weight update rule for complex signals is given by:

$$
\begin{align}
  y &amp;= \vec w_{n}^H x\\
  e_n &amp;= \left|y\right|^2 - 1 \\
  \vec w_{n+1} &amp;= \vec w_{n} - \mu \vec \nabla e_n \\
           &amp;= \vec w_{n} - \mu' e_n \Re\left\{y^\star x\right\}
\end{align}
$$

This can be implemented with:

for i = 1:N*M-Le,
    x = r(i:Le+i-1);
    y = w'*x;
    yd(i)=y;
    e(i) = abs(y)^2 - 1;
    w = w - mu * e(i) * real(conj(y) * x);
end


Plotting the results (showing the real part and imaginary parts in separate graphs) should give you something like the following: