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!


1 Answer 1


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.

enter image description here

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!

enter image description here

  • $\begingroup$ Thank you very much for this! It has been extremely helpful and I am currently trying to implement it on the microcontroller that I'm using. That approach to calculating the magnitude will surely save me some very valuable processing time. I will let you know how this turned out once I implement it. Thanks again! $\endgroup$
    – H A
    Mar 13, 2017 at 6:35
  • $\begingroup$ Super- Look into the Cordic Rotator too as that lends itself well to a microcontroller implementation for doing any trigomentric operations, including the complex rotation I first introduced (if you can't or don't want to multiply by $e^{j\omega t}$ directly, assuming you have the additional clock rate to allow for the iteration per cylce needed. en.wikipedia.org/wiki/CORDIC $\endgroup$ Mar 13, 2017 at 12:51
  • $\begingroup$ Thanks again! I'll also look into that. So far I have it all simulated on Matlab and I'm sure this will save me lots of processing time, compared with the FIR filters I've been using. $\endgroup$
    – H A
    Mar 14, 2017 at 7:10
  • $\begingroup$ Hi Dan thanks again for your help. I implemented the accumulate and dump filter and I've gotten great results. It requires around 5% of the operations I was doing for the FIR filters. As you mentioned, I ran several tests with an increasing number of ADC samples and the noise the noise decreases by sqrt(N). I have two questions that I hoped you could provide some insight with: 1) I extrapolated the results to obtain an expression Noise magnitude = a + b/sqrt(N). This tells me that a is the minimum amount of noise that my system will have, right? How can I interpret b? $\endgroup$
    – H A
    Mar 20, 2017 at 22:40
  • $\begingroup$ 2) How could I calculate the passband width of an accumulate and dump filter based on the amount of ADC samples? Thanks again! $\endgroup$
    – H A
    Mar 20, 2017 at 22:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.