In order to find the THD, we need to know the components that exist at each odd harmonic and the component that exists at the fundamental. We know our fundamental is at k=1 which means:
Since the rest of the sinusoidal terms in the square wave equation are harmonics of the fundamental, we know that:
By setting the initial value of ‘k’ back to one, we put the harmonic expression in a form with a known finite solution :
Plugging all of this into our THD equation above will yield a finite expression for the THD of a square wave:
Now, let’s say we want to multiply our square wave with a sine wave of the form . It doesn’t really matter what the amplitudes are here because they are multiplied by every single harmonic, so they are cancelled in the THD ratio. When we multiply with our defined square wave equation defined earlier, we will get the following:
That particular equation will go on to the infinith harmonic. Because of that, we can represent it as an infinite sum as follows:
Now that we have an expression, we can calculate the THD. Here, instead of our fundamental being at , it is actually at DC due to the demodulation. If we ignore the term (since it is in front of all of the harmonics, including DC), we see that our fundamental is at ‘1’ and all of the harmonics are contained within the infinite sum. Therfore, we can say our demodulated sine wave’s THD is equivalent to the coefficient of the infinite sum of . Now, lucky for us, said infinite sum also has a finite solution like the square wave did . Thus:
So the THD got worse! Well, does this actually make sense? Let’s look at it graphically:
I know the fourier transform of a sine wave is . For simplicity, instead of multiplying by a pure square-wave, let’s just say we only have the the harmonics ‘1’, ‘3’, and ‘5’. Given that, and normalizing the amplitudes to an arbitrary value of ‘A’, we can plot these signals in the frequency domain below.
Given this, if I convolve these two waveforms (remember, multiplication in the time-domain is convolution in the frequency-domain!), I will get the following plot (the math is very straight-forward so I’ll omit it here).
Now, if we sum up the squares of the harmonics and divide by the square of the component at DC, we will have . When we take the square root, we get . Which is larger than our exactly calculated value by despite only including two harmonics outside of the fundamental. As more harmonics are added into the square-wave approximation, the previous harmonic will decrease in amplitude (ie. the component at will go from to ).
I just thought this was an interesting result and decided to share it; hopefully someone, somewhere learned… something!
 I. V. Blagouchine and E. Moreau, “Analytic Method for the Computation of the Total Harmonic Distortion by the Cauchy Method of Residues” in IEEE Transactions on Communications vol 59, no. 9, pp. 2478-2491 Sept. 2011.  D. Zwillinger, Table of Integrals, Series, and Products, 8th ed. Oxford, UK: Elsevier, 2015, pp. 10-11