1
$\begingroup$

I am looking at FFT on events over the course of a year (~365 days). My data is the number of events in a day. I assign occurrence as a pulse with the amplitude reflecting the number of events. I see this as a pulse train because the events are essentially independent.

The following code includes functions that generate a signal as an array of amplitudes. I allow for varying periodicity of the pulses and the randomization of amplitudes. There is also a plotting function.

import numpy as np
import pandas as pd
from scipy.fft import fft, fftfreq, fftshift, ifft
import random
import holoviews as hv
hv.extension('bokeh')

def sig_spect(freq_df):

    ffts = hv.Scatter(freq_df.loc[freq_df.freq>0], kdims='freq', vdims=['spreal', 'period']
                      ).opts(color='red', xlabel='Freq', ylabel='Power',
                             default_tools=['hover', 'reset', 'box_zoom'], logx=False,
                             width=500, size=5)
    stems = hv.Spikes(freq_df.loc[freq_df.freq>0], kdims='freq',vdims=['spreal'])
    spectrum = hv.Overlay([stems, ffts])

    act_scat = hv.Scatter(signal_df).opts(width=300, default_tools=['hover', 'reset', 'box_zoom'])
    act_stem = hv.Spikes(signal_df)
    actual = hv.Overlay([act_scat, act_stem])

    return hv.Layout([actual, spectrum]).cols(2).opts(shared_axes=False)

def choice(i, spacing):    
    for space in spacing:
        if i%space == 0:
            return True
    return False

def sig_gen(num_samples = 365, spacing=[2], max_amp=1):

    # time in days
    t = np.arange(int(num_samples))
    # Amplitude and position of pulse. Amplitude here is 0 or 1 but can generate random values
    # Position here is every 7th day
    signal = [random.randint(1,max_amp) if (choice(i, spacing)) else 0 for i, x in enumerate(t)]
    
    # finds a multiple of K closest to and greater than a given number N
    def findNum(N, K):
        rem = (N + K) % K;

        if (rem == 0):
            return N
        else:
            return (N + K - rem)

    add_zero = findNum(num_samples, max(spacing)) - num_samples
    t = np.arange((num_samples + add_zero))
    signal = np.pad(signal,(0,add_zero))
    return t, signal

I run the following code with 1 event every 7 days:

t, signal = sig_gen(spacing=[7])

sp = fft(signal)
freq = fftfreq(t.shape[-1])

freq_df = pd.DataFrame({
    'freq': freq,
    'spreal': np.abs(sp),
    'period': 1 / freq
})
signal_df = pd.DataFrame({
    't' : t,
    's' : signal
})
sig_spect(freq_df)

It results in a spectrum I expect. enter image description here

If I take the same period but randomly vary the amplitude between 0 and 10 it is still what I expect: enter image description here

If however I combine two periods, 4 days and 7 days the results are less clear:

t, signal = sig_gen(spacing=[4, 7], max_amp=1)

enter image description here

The spectrum is even noisier if I introduce varying amplitudes and additional periods. My question is:

  • if I know what periods to expect then I can easily identify their components in the spectrum. If I don't know what to expect then how would I do so? Note that in this case I can usually assume a non-integer period is not possible but what happens in a more general case.

Note 1: I'm currently padding for the largest period. That may introduce a problem if it isn't suitable for all the periods. How would you pad to accommodate all periods of interest?

$\endgroup$
4
  • $\begingroup$ This seems like PRI estimation ([multi-]pulse repetition interval estimation), a problem often encountered in radar systems. $\endgroup$
    – Peter K.
    Dec 16, 2021 at 18:36
  • 1
    $\begingroup$ @PeterK. I looked that up and I guess it does look like that. $\endgroup$
    – MikeB2019x
    Dec 16, 2021 at 20:14
  • $\begingroup$ I didn't review your code but my intuition would be zero insert in time such that you have plenty of zero samples in between your events, which is basically sampling at a higher rate. This combined with additional zero padding to the degree you want to interpolate the frequency result should give satisfactory results, to the extent that the "noise" in the variation from sample to sample does not exceed the signal level itself. $\endgroup$ Dec 17, 2021 at 1:01
  • $\begingroup$ I have tried increasing padding (zeros) just at the end and at both ends of the signal. The addition changes peak heights a small amount but the lower period component always dominates (naturally) and if I weren't aware of what I was looking for I would likely miss the higher period component. $\endgroup$
    – MikeB2019x
    Dec 17, 2021 at 13:53

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.