They say a picture is worth a thousand words. And a handful of animated pictures...
All of these from Wikipedia. You can see these in context in the articles for Fourier transform, square wave, and Fourier series.
As it turns out, any periodic waveform can be represented by a Fourier series, although to represent it exactly may require an infinitely long series. In engineering practice, however, above a certain point continuing the series adds no additional useful precision, and we can call a finite series "good enough".
Direct answers to your questions:
How is it possible for something to generate a waveform across all frequencies by just looking at a waveform from one frequency?
I think your mental hangup here is "one frequency". The point of Fourier analysis is that a square wave (say, a 60 Hz square wave) is not just one frequency. I mean it is, in the sense that it's just a 60 Hz square wave. But we can also decompose that square wave into a number of sinusoidal waves, each with their own amplitude and phase.
Also, if harmonics are the changes in gain for other frequencies, why do I care about them if I'm running at a single constant frequency?
Because, through the lens of Fourier analysis, you aren't actually at a single constant frequency unless your waveform is a completely undistorted sine wave.
How can these changes from other frequencies affect something from another frequency?
In a linear and time invariant system, they can't. This is formally described in LTI system theory. Each frequency component passes through the system completely independently. We can consider what happens for one frequency, then separately consider what happens for another frequency. Adding the two together, we get the same result as if the two frequencies passed through the system simultaneously. This makes Fourier analysis a very useful analytical tool.
Also, is an RLC band-stop filter a good way of cleaning up the waveform by simply blocking out the frequency range where the harmonics are? How would you block out multiple frequency ranges since there's more than one bad harmonic?
If the problematic harmonics are in just one band, yes. In a power system, it's likely that you will be dealing with 60 or 50 Hz with some degree of harmonic distortion, meaning it will also contain frequency components at integer multiples of the fundamental. If you can design a filter that allows that fundamental frequency to pass, but attenuates all higher harmonics, then no matter how distorted the input, the output will be a pure, undistorted sine wave. In a power system you might also care about the phase of voltage and current (power factor), and there are filters to fix that too.
You're confusing number of samples with resolution.
The 16-bit transform means each sample is 16-bits. In other words, the samples can be -32768 to 32767. The number of samples (N_WAVE) taken is completely irrelevant to the 'bitness' of the FFT. The 16-bits is referring to sample size, not the number of samples being taken. So the 1024 sample count and the resolution of those said samples have no relationship to each other, and 'where it came from' is, well, from the number of samples you wanted to take of a wave. This is determined by the nyquist frequency. You can just as readily use an N_WAVE of 256 if you wanted.
As for why your code isn't working, you're FFT is probably overflowing like crazy. The PIC doesn't magically change multiplication of 2 16-bit values into a 32 bit value like it does for turning 2 8-bit values into a 16-bit value. You need more than 16-bits of storage if you're multiplying 16-bit numbers. Even if your dynamic range is just 12 bits, you're scaling it, and if you didn't, well, anything except 0 multiplied by some of the values in your lookup table (which span the entire signed short value range) will overflow.
In the fix_fft function, use ints instead of shorts for the two arrays and m. This will force 32 bit multiplication and will resolve any overflow issues.
Honestly, this is probably not a question suited to electronics, so you might get a better answer in a programming related stackexchange. But I think this will help. Good luck!
Best Answer
√(I2+Q2) just nets you the absolute value (amplitude) of a sample. Doing that for every sample gives you the envelope (red trace) of the signal
which isn't useful for what you want to do. Certainly not for computing distances or velocities.
Instead, treat the I, Q samples as complex numbers of the form I + Qj (where j is the imaginary unit) and compute the discrete fourier transform ("FFT") of that complex-valued data. Only then compute the absolute value (√(I2+Q2)) of each bin of the resulting spectrum.
Each spike in the spectrum then (ideally) corresponds to a target. With fixed targets, the frequencies of the spikes are proportional to the range to the respective targets while the amplitudes reflect the strength of the echoes.
With moving targets, the frequency of a spike is dependant on a combination of the radial velocity and range. With FMCW radar you are always measuring velocity as well, even when you'd rather not. This is because the doppler shift interferes with the distance-induced frequency difference being measured.
The frequency offset due to the doppler shift is unaffected by the direction and rate of a transmitted chirp (whether the frequency increases or decreases with time and how fast), while the range-induced frequency offset very much is affected. In the extreme case of transmitting a constant frequency, you are only sensing the velocity (via the doppler shift), hence such a radar is termed doppler radar.
This can be exploited to separate the range from the velocity by constantly changing chirp directions and rates and then performing some fancy math. How to do this in practice is outside the scope of an answer (and my brain).
Assuming that you are using floating point math and that you don't have a significant DC offset, I'd just map the I and Q samples from the input range [0, 216) to [-1, 1]. In essence, I = I/215-1, Q = Q/215-1.