fir - Obtaining phase and magnitude of 2 spectral components after ADC

I am currently implementing a digital filter to obtain the output of a PID controller. The control loop is carried out by applying a 1kHz sine wave at the input of the system and reading the output signal. The signal is distorted by the system and the output signal contains a 2nd harmonic component of the 1kHz input sine wave. I am currently applying FIR filters to obtain the phase and magnitude of the fundamental and harmonic tones but in order to reduce the noise I have to use very high order FIRs (around 1000 taps) and it is a very slow process. I would like to know if a matched filter would be a better approach, or if there are any other suggestions. Thanks!

Answer

To obtain the phase and magnitude of a known frequency, I recommend multiplying each sample by $exp^{-j\omega t}$ where $\omega$ is $2\pi 1000$ Hz for the case of the fundamental (and $2\pi n 1000$ for each harmonic n). The output can then be "accumulated and dumped" over a number of samples depending on how much noise is present, traded with possible offset errors in the actual frequency, and how stationary the process is. For every N samples accumulated, the noise if white will go down relative to the accumulated output by $\sqrt{N}$. If this measurement is being used in the control loop itself then there is the additional complexity of ensuring that the accumulation time is significantly shorter than your loop response time or this process will become part of your loop filter, which you probably don't want! (In this case I would likely let the loop filter take care of all the filtering.)

This "accumulate and dump" is equivalent to a moving average filter followed by a decimator, so a 1 ms accumulation time would be very attractive as the nulls would all be spaced by 1 KHz in that case (conveniently nulling out all the harmonics for all cases)! Further this "correlator" which is what this is would have an amplitude vs frequency offset sensitivity given by the Sinc frequency response of the moving average filter.

Then with the accumulate and dump output (I and Q using the complex method I described) phase, frequency and magnitude estimates can be easily obtained with minimum processing which can additionally be further averaged to increase measurement SNR (if the system is sufficiently stationary), or include in a correction loop which I suspect the OP may be trying to do.

Here is a common magnitude algorithms for efficiently estimating magnitude from I and Q samples (if one doesn't want to or can't calculate $\sqrt{I^2+Q^2}$ directly): http://dspguru.com/dsp/tricks/magnitude-estimator/.

The phase information is the phase of the I,Q sample since we have translated the signal of interest to baseband with a complex rotation (so in that case we need to compute $ATAN2(Q,I)$ which can be done with a look-up table or the Cordic Rotator if an iterative approach is allowed.

NOTE: If the intention was to implement a phase lock loop, then this becomes even simpler in that a multiplier can be used as the phase detector (as the output of a multiplier would be proportional to the cosine of the phase between the two input signals), knowing that the loop will drive the phase detector to it's linear region at the 90° crossing, and in this case a real 1KHz signal can be multiplied by the test signal.

For efficient frequency estimation please see the figure below that shows how a frequency estimate can be extracted from this simple operation on two consecutive complex samples:

$$Freq \varpropto I[n]Q[n-1]-I[n-1]Q[n]$$

This comes about because the complex conjugate multiplication of any two vectors will be a vector with a phase that is equal to the phase between the two vectors (so for small angles and normalized vectors, the Q of that vector is the phase!). This is known as a cross-product discriminator for estimating phase (and in the case when the two vectors are the same signal delayed in time, this would be change in phase over a change in time which is frequency). Any such frequency measured would be the frequency offset from the 1KHz, or whichever harmonic was selected.

This plot below shows how the phase is derived (and here the vectors $V_1$ and $V_2$ do not relate to the same vectors in the previous plot). The dashed line shows the conjugate of $V_1$. When two vectors are multiplied, their magnitudes multiply and their angles add, resulting in the product vector drawn in red that I will call $V_3$. If $V_1$ and $V_2$ were normalized to be 1, and if the angle between $V_1$ and $V_2$ was small, then the imaginary portion of $V_3$ (Q) would be equal to the angle with no arctan calculation necessary!