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.
Instantaneous
to mean "only dependent on the current time domain sample". $\endgroup$