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

[mantzaris]

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?

[ssfrr]

This isn’t a complete answer to your question, but maybe some helpful pointers.

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, φ
        # assumes ω is in radians/sample and φ is in radians
        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.

[ssfrr]

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.

[mantzaris]

@ssfrr

Once you have the data you’re looking for (say freqs , amps , phases are each a Vector of respective values)

How do I obtain those 3 values?

[ssfrr]

@mantzaris

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.).

[mantzaris]

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

[StevenSiew]

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.

almost forgot about ModelFunc

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