I am trying to calculate the instantaneous power of a 1-d signal within a particular frequency band in an on-line, real-time, streaming application. I have tried using sliding window FFTs but these require computing an average power over a window rather than the instantaneous power at each incoming sample. Wavelet transforms also appear to not be suitable for on-line or real-time instantaneous power because of edge effects in the analysis window. I have considered using Hilbert-Huang Transforms, but these seem to have the same issue as wavelet transforms with edge effects and compute time. RMS analysis is fastest, but still returns average power over a window of samples rather than the power at a particular sample.
Problem statement: Given a data stream
data_stream = x_0, x_1, x_2, x_3..., x_n where
x_i is a single, real-valued number, how to compute the the instantaneous power of the signal at
x_i for a particular frequency band?
Attempts at FFT and RMS solutions (pseudocode):
- Sliding window FFT:
window_length = 50 # number of samples required to capture 2 cycles of lowest frequency in frequency band of interest low_f = 0.25 # 0.25 Hz lowest frequency of interest high_f = 0.45 # 0.45 Hz highest frequency of interest band = [low_f, high_f] # band of interest includes all components of signal from 0.25-0.45 Hz start_index = 0 end_index = start_index + window_length while end_index < len(data_stream): window = data_stream[start_index: end_index] fft = fft(window) # compute fft on window power = fft[band] start_index += 1 end_index += 1
Here, I keep the window length constant at 50 samples, and move the window one sample forward each time I receive a new sample. The band
power is the average power of the band from all 50 samples in the window
data_stream[start_index: end_index] rather than the instantaneous power of the signal at
data_stream[end_index - 1].
- Sliding window RMS:
filter = bandpass_filter(low_f, high_f) window_length = 50 # number of samples required to capture 2 cycles of lowest frequency in frequency band of interest start_index = 0 end_index = start_index + window_length while end_index < len(data_stream): window = data_stream[start_index: end_index] filtered_window = filter(window) # reduce signal to only frequencies of interest power = RMS(window) start_index += 1 end_index += 1
Here, I bandpass filter each window of the
data_stream to keep only the frequency band of interest. I keep the filtered window length constant at 50 samples, and move the filtered window one sample forward each time I receive a new sample. The band
power is the average power of the band from all 50 samples in the filtered window
filtered(data_stream[start_index: end_index]) rather than the instantaneous power of the signal at
data_stream[end_index - 1]. This sliding window RMS approach should give me approximately the same results as the sliding window FFT as FFT power is ~ RMS.
- "Instantaneous Power":
online_filter = online_bandpass_filter(low_f, high_f) # see https://stackoverflow.com/questions/40483518/how-to-real-time-filter-with-scipy-and-lfilter current_index = 0 while current_index < len(data_stream): current_sample = data_stream(current_index) # get new sample filtered_sample = online_filter(current_sample) power = filtered_sample ** 2 current_index += 1
Here I use an on-line bandpass filter (as discussed in this post)on each sample as it is received. Assuming that
filtered_sample represents only the frequency band of interest, I then square the value of that sample to get the "instantaneous power".
As far as I understand, the way real-time oscilloscopes / spectroscopes work is computing the FFTs on a sliding window of incoming samples, and then presenting the FFT PSD for each window as the "real-time" PSD.
For example, if a window of
2.0s is required to capture all frequencies of interest, and the sampling rate is
2 Hz, then the oscilloscope /spectroscope computes the FFT PSD on the samples from
[t = -2.0s, t = 0.0s] and presents that PSD as the "instant, real-time" PSD at time
t = 0.0s. Another sample is obtained, and the window shifts to
[t = -1.5s, t = 0.5s]. FFT PSD is computed on the new window, and that PSD is reported as the "instant, real-time" PSD at time
t = 0.5s.
In reality, the PSD computed from
[t = -2.0s, t = 0.0s] is the PSD over the interval rather than the PSD at
t = 0.0s, and likewise for the PSD computed at the next sample acquisition from
[t = -1.5s, t = 0.5s].
Am I understanding this correctly? Is there a better way to calculate the power of a signal at the latest available data point while new data arrives in a on-line / streaming manner? My attempt #3 ("Instantaneous Power") seems to be the closest, but I'm not sure if it's theoretically sound. The output is super high variance, noisy, and when averaged, not equivalent to the FFT output.
Thanks for your thoughts!
EDIT: Edited wrt to comments and previous answer.