Robin Popelka - 3 months ago 7

Python Question

I am just fighting with modulation of sinus wave.

I have got a frequency (from messured data - changing in time) and now I need to plot sinus wave with coresponding frequency.

The blue line are just plotted points of real data and the green is what I did till now, but it does not corespond with real data at all.

The code to plot sin wave is bottom:

`def plotmodulsin():`

n = 530

f1, f2 = 16, 50 # frequency

t = linspace(6.94,8.2,530)

dt = t[1] - t[0] # needed for integration

print t[1]

print t[0]

f_inst = logspace(log10(f1), log10(f2), n)

phi = 2 * pi * cumsum(f_inst) * dt # integrate to get phase

pylab.plot(t, 5*sin(phi))

Amplitude vector:

[2.64, -2.64, 6.14, -6.14, 9.56, -9.56, 12.57, -12.57, 15.55, -15.55, 18.04, -18.04, 21.17, -21.17, 23.34, -23.34, 25.86, -25.86, 28.03, -28.03, 30.49, -30.49, 33.28, -33.28, 35.36, -35.36, 36.47, -36.47, 38.86, -38.86, 41.49, -41.49, 42.91, -42.91, 44.41, -44.41, 45.98, -45.98, 47.63, -47.63, 47.63, -47.63, 51.23, -51.23, 51.23, -51.23, 53.18, -53.18, 55.24, -55.24, 55.24, -55.24, 55.24, -55.24, 57.43, -57.43, 57.43, -57.43, 59.75, -59.75, 59.75, -59.75, 59.75, -59.75, 59.75, -59.75, 62.22, -62.22, 59.75, -59.75, 62.22, -62.22, 59.75, -59.75, 62.22, -62.22, 62.22, -62.22, 59.75, -59.75, 62.22, -62.22, 62.22, -62.22, 59.75, -59.75, 62.22, -62.22, 62.22, -62.22, 62.22, -62.22, 59.75, -59.75, 62.22, -62.22, 59.75, -59.75, 62.22, -62.22, 59.75, -59.75, 59.75]

Time vector for real data:

[6.954, 6.985, 7.016, 7.041, 7.066, 7.088, 7.11, 7.13, 7.149, 7.167, 7.186, 7.202, 7.219, 7.235, 7.251, 7.266, 7.282, 7.296, 7.311, 7.325, 7.339, 7.352, 7.366, 7.379, 7.392, 7.404, 7.417, 7.43, 7.442, 7.454, 7.466, 7.478, 7.49, 7.501, 7.513, 7.524, 7.536, 7.547, 7.558, 7.569, 7.58, 7.591, 7.602, 7.613, 7.624, 7.634, 7.645, 7.655, 7.666, 7.676, 7.686, 7.697, 7.707, 7.717, 7.728, 7.738, 7.748, 7.758, 7.768, 7.778, 7.788, 7.798, 7.808, 7.818, 7.828, 7.838, 7.848, 7.858, 7.868, 7.877, 7.887, 7.897, 7.907, 7.917, 7.927, 7.937, 7.946, 7.956, 7.966, 7.976, 7.986, 7.996, 8.006, 8.016, 8.026, 8.035, 8.045, 8.055, 8.065, 8.075, 8.084, 8.094, 8.104, 8.114, 8.124, 8.134, 8.144, 8.154, 8.164, 8.174, 8.184, 8.194, 8.20]

So I need generate sinus with constant amplitude and following frequency:

[10.5, 16.03, 20.0, 22.94, 25.51, 27.47, 29.76, 31.25, 32.89, 34.25, 35.71, 37.31, 38.46, 39.06, 40.32, 41.67, 42.37, 43.1, 43.86, 44.64, 44.64, 46.3, 46.3, 47.17, 48.08, 48.08, 48.08, 49.02, 49.02, 50.0, 50.0, 50.0, 50.0]

Answer

You can try to match you function with something sine- or actually cosine-like, by extracting estimates for the frequency and the amplitude from your data. If I understood you correctly, your data are maximums and minimums and you want to have a trigonometric function that resembles that. If your data is saved in two arrays `time`

and `value`

, amplitude estimates are simply given by `np.abs(value)`

. Frequencies are given as the inverse of two times the time difference between a maximum and a minimum. `freq = 0.5/(time[1:]-time[:-1])`

gives you frequency estimates for the mid points of each time interval. The corresponding times are thus given as `freqTimes = (time[1:]+time[:-1])/2.`

.

To get a smoother curve, you can now interpolate those amplitude and frequency values to get estimates for the values in between. A very simple way to do this is by use of `np.interp`

, which will do a simple linear interpolation. You will have to specify at which points in time to interpolate. We will construct an array for that and then interpolate by:

```
n = 10000
timesToInterpolate = np.linspace(time[0], time[-1], n, endpoint=True)
freqInterpolated = np.interp(timesToInterpolate, freqTimes, freq)
amplInterpolated = np.interp(timesToInterpolate, time, np.abs(value))
```

Now you can do the integration, that you already had in your example by doing:

```
phi = (2*np.pi*np.cumsum(freqInterpolated)
*(timesToInterpolate[1]-timesToInterpolate[0]))
```

And now you can plot. So putting it all together gives you:

```
import numpy as np
import matplotlib.pyplot as plt
time = np.array([6.954, 6.985, 7.016, 7.041, 7.066, 7.088, 7.11, 7.13]) #...
value = np.array([2.64, -2.64, 6.14, -6.14, 9.56, -9.56, 12.57, -12.57]) #...
freq = 0.5/(time[1:]-time[:-1])
freqTimes = (time[1:]+time[:-1])/2.
n = 10000
timesToInterpolate = np.linspace(time[0], time[-1], n, endpoint=True)
freqInterpolated = np.interp(timesToInterpolate, freqTimes, freq)
amplInterpolated = np.interp(timesToInterpolate, time, np.abs(value))
phi = (2*np.pi*np.cumsum(freqInterpolated)
*(timesToInterpolate[1]-timesToInterpolate[0]))
plt.plot(time, value)
plt.plot(timesToInterpolate, amplInterpolated*np.cos(phi)) #or np.sin(phi+np.pi/2)
plt.show()
```

The result looks like this (if you include the full arrays):