I have Analog to Digital converter from which I get data of measurements of 1kg load. The aim is to model the simple spring-damper system, so that after entering correct initial data (x[0] = a, x'[0] = b) the model will give almost the same data as the original data from converter.

It is my first modeling attempt of a such real system. The problem is that the data from converter contains the 50Hz noise (electricity frequency) and if I use FIR filter I am afraid that the initial values needed for modeling the same wave data will be corrupted. How this process is done usually? As for modeling equation I use mx'' = -kx - c x'. And I have made numerical solution for this equation.

Edit 1: How about not filtering but modeling the system with taking into account the electrical noise?

  • $\begingroup$ What input are you driving the system with? $\endgroup$
    – Jason R
    Feb 11 '12 at 1:19
  • $\begingroup$ @Jason R: Input is from A/D converter, but with some convertion to the displacement valuse of x[t] (not raw input of A/D converter) $\endgroup$
    – maximus
    Feb 13 '12 at 6:48
  • $\begingroup$ No, I mean what input signal are you forcing the load with? $\endgroup$
    – Jason R
    Feb 13 '12 at 13:56
  • $\begingroup$ @JasonR Sorry, still can't understand your question. What input signal are you talking about? I have signal from scales, thats all signals I have. $\endgroup$
    – maximus
    Feb 14 '12 at 2:30
  • $\begingroup$ How does the simulation start? Do you apply a force to the mass? You have to do something to the system to get it to start moving at all. If you can control the force that you put in, then you may be able to obtain an output displacement signal that is more robust to filtering around the 50 Hz electrical signal frequency. $\endgroup$
    – Jason R
    Feb 14 '12 at 13:29

The problem is that just taking the x[0] and x'[0] of your observed data and plugging them into the differential equation will be extremely sensitive to any kind of noise ; not only the 50Hz hum. You are trusting those two samples (rather than the entirety of your observed data) for your estimation, and you will get incorrect results from that, whether the 50Hz noise is present or not.

How to do it right? First, get rid of the 50Hz noise with a notch filter tuned to 50Hz and maybe to a few harmonics.

Then, what you want to do is fit the parameters $A$ and $\phi$ of a $x(n) = A e^{-\delta n} sin(\omega n + \phi)$ model to your data $y(n) = x(n) + noise$. You already know $\delta$ and $\omega$ since they depend on your known system parameter (damping and oscillation frequency). What you really want to know is the best estimate of $A$ and $\phi$ to reconstruct the "theoretical" part (in proper speak: signal subspace projection) of your observation.

Here is how you can estimate $A$ and $\phi$ without lame methods involving looking at peaks and zerocrossing:

You can write your signal model as $x(n) = Re (\alpha z^n) = \frac{1}{2}\alpha z^n - \frac{1}{2} \bar{\alpha} \bar{z}^n$ where $z$ is the known complex pole $z = e^{-\delta + j\omega}$, $\alpha$ is the complex amplitude $\alpha = Ae^{j \phi}$ and $\bar{z}$ denotes the conjugate.

Lay out a $N\times2$ matrix where N is your number of observation, the first column being $\frac{1}{2} [1, z, \ldots, z^N]$ and the second $-\frac{1}{2} [1, \bar{z}, \ldots, \bar{z}^N]$. Do a least square regression of that matrix (which is a basis of your signal subspace) onto the data matrix $[y(1), \ldots, y(N)]$. This will give you LS estimates of $\alpha$ and $\bar{\alpha}$ as the two solutions. Take the modulus and angle of the first to get $A$ and $\phi$. Plug this back into $A e^{-\delta n} sin(\omega n + \phi)$ to get the best fit theoretical model from your data.

# Estimate amplitude and phase of an exponentially damped model
# given observations y and system parameters delta and omega
n = numpy.arange(y.shape[0])
z = numpy.exp(-n * delta + n * omega * 1j)
z_conj = numpy.vstack((0.5 * z, -0.5 * numpy.conj(z))).T
s, _, _, _ = numpy.linalg.lstsq(z_conj, y)
print "Estimated amplitude: ", numpy.abs(s[0])
print "Estimated phase: ", numpy.angle(s[0])

If you want to also estimate $\delta$ and $\omega$ the procedure is a tad more complex and involves an EVD or SVD of the correlation matrix of your data to get a signal subspace basis, one more EVD to get the estimate of $z$, then the procedure sketched above to retrieve $A$ and $\phi$. You will have a robust estimate of all the parameters of the theoretical model from your data. Search for ESPRIT or MUSIC.

Note that this approach can be used to estimate the phase / frequency of a mixture of any number of damped sine waves. So for example, you could use this to retrieve the amplitude and phase of the 50Hz hum itself. Just put two additional columns in the matrix used for the least square procedure, with powers of $e^{j\omega_{noise}}$ and their conjugates. This will give the best fit amplitudes and phases of an exponentially damped sine wave (described by omega and delta) + 50 Hz hum model.

  • $\begingroup$ I have posted an example of estimation of the 50Hz hum amplitude + model parameters: gist.github.com/1788715 $\endgroup$ Feb 10 '12 at 10:48
  • $\begingroup$ For example in case if I plase a 3kg weight, then the frequency of the system will be very close to the electrical noise frequency, and what I am afraid of is that the result after notch filtering will be not adequate for modeling, because I filter out not only the noise. $\endgroup$
    – maximus
    Feb 13 '12 at 6:50
  • $\begingroup$ By the way very nice solution using linear LSM instead of non-linear, can you provide any link or book you used to learn this methods? $\endgroup$
    – maximus
    Sep 11 '12 at 4:45

If the interference from line power isn't so large that it is distorting your ADC's output (i.e. the large 50 Hz component isn't driving the dynamic range of the overall sampled signal), then you could design a notch filter to remove the narrowband line frequency component. Choosing the best way to implement such a filter depends upon the computational resources that you have, the type of arithmetic that you will use, other requirements on the filter (how wide can the notch be? how flat does the passband have to be?), etc. Here is a description of how to design a second-order IIR notch filter.

If the frequency drifts over time sufficiently enough that you can't just place a static notch in the spectrum that is wide enough to cover its entire range, then you may need to consider an adaptive approach. This is significantly more complicated, but it is possible to build a periodic interference suppressor using a least-mean squares (LMS) filter.

  • $\begingroup$ I have checked, in my case the problem is that the noise frequency and the real data frequency are very close to each other, how to deal with it regarding your answer? $\endgroup$
    – maximus
    Feb 8 '12 at 2:08

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.