Electronic – Frequency of signal

You are confusing normalized frequency with real frequency (Hz). Wikipedia Normalized Frequency.

The equation for the fundamental frequency (\$f_0\$) given number of delay units in the "string" (\$M\$) should be \$f_0=\frac{F_s}{M}\$, as the website states. Note that the Karplus Strong algorithm (formulated for a plucked string) is no more than a feedback comb filter initialized to a noise burst, which has frequency response: \$H(\omega)=\frac{1}{1-\alpha e^{-j \omega M }}\$, where \$\alpha\$ is the scaling factor from your original difference equation.

Here is the magnitude plot with \$\alpha=.5\$ and \$M = 10\$:

Wolfram Alpha Plot

Note that the peaks are at \$\omega = 0, \pi/5, 2\pi/5...\$ This implies the fundamental frequency is at \$\omega_0 = \pi/5\$, which corresponds to the "real" frequency of \$44100/10 = 4410\$ Hz (assuming 44100 samples / sec). Remember that the fifth harmonic in this case is \$\omega_5 =5\pi/5 = \pi\$. This corresponds to what we expect, the fifth harmonic is 22050 Hz, the highest frequency we can represent with a sample rate of 44100 samples / sec.

Think about what happens in the frequency domain: over time, the initialized noise burse (which looks fairly uniform in frequency) gets filtered by the peaked harmonic response of the filter, which produces the harmonics you expect for a plucked string.

In general, Fourier coefficients are given by the formulas:

$$a_n = \frac{2}{T}\int_{x_0}^{x_0+T}g(x)\cos\left(\frac{2\pi nx}{T}\right)dx\\ b_n = \frac{2}{T}\int_{x_0}^{x_0+T}g(x)\sin\left(\frac{2\pi nx}{T}\right)dx$$

For any integrable \$g(x)\$ on the interval \$[x_0, x_0+T]\$. For periodic functions this interval can correspond to the period of the function, such that the Fourier series will be perfectly equal to the function itself. But the case of your example is very simple - it's just a sum of three cosines, so the \$a_n\$ terms can be calculated directly and there would be three of them. \$b_n\$ would equal to zero, as there is no sine terms.