# How to obtain the frequency/phase/power of a signal and then produce a waveform of arbitrary length from those parameters?

I have a waveform signal `X1` in the time domain that is long enough for all of the frequency components to have completed multiple periods, and want to derive the `frequency`, `phase` and `amplitude/power` components.

With those components derived, I would like to then subsequently construct the time domain waveform for an arbitrary selected length using those components. How can I do this in Julia?

Can I print/display the list of the freq/phase/amplitude components of the signal to understand its composition to be interpreted easily?

The first thing Iâd do is take a look at the spectrum with `plot(abs.(rfft(X1)))` If the components are sufficiently far apart from each other you should seem them clearly as peaks. Note though that in general the components wonât lie in the center of the FFT bins, so youâd want to do some kind of interpolation (e.g. cubic) to better estimate the actual peak.

Do you expect the components to be harmonics of an underlying fundamental? If so then that fundamental gives you the period where the whole signal repeats. You can use `xcorr(X1, X1)` to look for the peak (though if the period is a substantial fraction of the signal length youâll run into problems near the edges because the cross-correlation rolls off). Once you have the period, you can truncate to that length and the components will lie in the middle of the FFT bins, so you donât need to worry about interpolation.

Once you have the data youâre looking for (say `freqs`, `amps`, `phases` are each a `Vector` of respective values), then creating the time-domain signal of length `L` is something like:

``````sig = map(0:L-1) do n
mapreduce(+, freqs, amps, phases) do Ď, a, Ď
a * cos(2Ď*Ď*n + Ď)
end
end
``````

Maybe a little over-clever with `map` and `mapreduce`, but this doesnât need to make any temporary arrays so should be pretty fast.

2 Likes

Actually another way to do the frequency estimation would be to use âinstantaneous frequencyâ (IF), which is usually defined as the phase derivative in each FFT bin. You can do a 1-sample finite-difference approximation by just taking the FFT of your signal as well as a 1-sample delayed version of the signal. For each FFT bin the phase difference between the two signals should be the frequency in radians/sample. Youâll probably want to apply a window function (e.g. `hanning` or `gaussian` from DSP.jl) to your signal to contain spectral leakage.

Thereâs a fancier version of the IF thing where instead of delaying by one sample for the second signal you multiply it by the derivative window, but Iâd need to refresh myself on the math - there may be some normalization you have to add.

I think both of these may be sensitive to noise, so you might want to watch out if your signal-to-noise is bad.

3 Likes

How do I obtain those 3 values?

Thatâs what the rest of the post was referring to. Youâll probably have to try some things to see what works best for your signal, depending on some of the issues I brought up above (signal-to-noise, signal length, how closely-spaced the components are, whether they are harmonics, etc.).

3 Likes

for example:
`x = sin.(2 .* (0:0.1:10*pi)) + sin.(0:0.1:10*pi)`
with `plot(abs.(rfft(x)))` produces;

I can see the centers for the 2 frequencies, but the âamplitudeâ and phi? Would the Hilbert work on a signal with multiple frequencies would the different modes not require separation first? From this plot the 3 variables of omega, a, and phi will be inferred ?

edit: `plot(fftfreq(x),abs.rfft(x))` (or `plot(fftfreq(x),abs.rfft(x/N))`) will give the magnitudes along the frequency bins and then using the function `angle(x)` at each frequency position where there is a large peak produces the phase. This delivers what is needed for subsequent reconstruction

This is everything I know about how to extract the signal from the time domain data. I am sorry that this is in Mathematica but the maths is all there.

``````ModelFunc[data_, mag_, freq_, angle_] := Module[{d},
d = data -
Table[
N[ mag * Sin[freq*t + angle], 14]    , {t, 0,
Length[data] - 1}];
Abs[X[d, freq]]
]
``````
``````data = {0, 0, 17, 4, 2, 2, 3, 16, 5, 6, 6, 7, 5, 4, 2, 2, 2, 4, 6, 7,
7, 7, 5, 4, 2, 1, 1, 2, 2, 3, 2, 2, 2, 2, 1, 2, 2, 2, 3, 2, 2, 1,
0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 2, 3, 3, 3, 4, 4, 4, 7, 10, 12,
10, 8, 6, 6, 7, 8, 10, 24, 36, 24, 10, 7, 6, 6, 7, 7, 8, 7, 8, 7,
8, 9, 9, 9, 9, 9, 11, 11, 13, 15, 17, 20, 21, 22, 22, 22, 21, 22,
24, 21, 18, 13, 10, 10, 11, 13, 13, 15, 17, 18, 25, 31, 28, 22, 23,
33, 38, 37, 31, 29, 30, 32, 34, 37, 46, 83, 81, 42, 33, 32, 36,
32, 23, 19, 27, 55, 72, 81, 71, 39, 31, 43, 94, 139};

PlotAndModelFreq[data, 64]
``````

You get something like this