Least square regression of complex amplitudes onto a synthetic signal subspace for an EDS+50 Hz hum model

gistfile1.py

import numpy
import pylab

delta = 0.01
omega = 0.05
omega_noise = 0.2
A = 5
A_noise = 0.2
phi = 0.8
N = 1000
n = numpy.arange(N)

# "Theoretical" part
x = A * numpy.exp(-n * delta) * numpy.cos(n * omega + phi)
y = x + A_noise * numpy.cos(n * omega_noise)

n = numpy.arange(y.shape[0])
z = numpy.exp(-n * delta + n * omega * 1j)
z_noise = numpy.exp(n * omega_noise * 1j)

z_signal = numpy.vstack((
    0.5 * z,
    -0.5 * numpy.conj(z),
    0.5 * z_noise,
    -0.5 * numpy.conj(z_noise))).T

s, _, _, _ = numpy.linalg.lstsq(z_signal, y)
print s
print "Estimated amplitude: ", numpy.abs(s[0])
print "Estimated phase: ", numpy.angle(s[0])
print "Estimated noise amplitude: ", numpy.abs(s[2])
print "Estimated noise phase: ", numpy.angle(s[2])

pylab.plot(y)
pylab.plot(numpy.dot(z_signal[:, :2], s[:2]).real)
pylab.savefig('figure.pdf')