Strange spectrum results from sim_knee?

[SusanL82]

Hi all,

I've been trying to use the neurodsp signal simulation functions to try to understand how aperiodic parameters affect some other spectrum measures from our in vivo experiments (changes in a power-ratio-based state space). I wanted to blame changes in knee values for some results, so I've used sim_knee with a variety of knee values and the same two slopes and with a sampling rate matching my real signals (6kHz).

I noticed that the power spectra of the simulated signals look very strange. I used a sliding window FFT for 2-second windows (50% overlap, for 60min of simulated data) and then averaged the spectrum for all windows. This is the method we use for the real data.

When I used somewhat randomly chosen values (exp1 = -0.5, exp2 = -1) and got this:

I then tried settings for knees and exponents that are much more in the range of what I see in my recordings (exp1 = -6, exp2 = -8), but I get the same odd repeating pattern with huge unexpected peaks:

Am I just doing something wrong with sim_knee? Or does the function only work with specific ranges of inputs?

[SusanL82]

The scripts I used for simulating my signals looks like this. I'm making a lot of small sections that I merge together in a later step. If I try to simulate an hour in one go, it first is extremely slow and then tends to crash at some point.

`# -*- coding: utf-8 -*-

# Set some general settings, to be used across all simulations
fs = 6000
n_seconds = 60
knee_list =[200,500,1000,1500,2000,5000,10000] #[25,50,100,200,500,1000,1500]
OutFolder = "D:/SimSigs"

# Import sim functions
from neurodsp.sim import (sim_knee)
from neurodsp.utils import set_random_seed
from scipy.io import savemat
import numpy as np

# Import utilities for plotting data
from neurodsp.utils import create_times

# Create a times vector for the simulations
times = create_times(n_seconds, fs)

for k in range(np.size(knee_list)):
    for l in range(1 ,60): 

        # Set the random seed, knee is only thing different for signals
        set_random_seed(l)

        # Simulate a knee signal, with specified exponents & knee
        exp1 = -6
        exp2 = -8
        knee_ap = sim_knee(n_seconds, fs, exponent1=exp1, exponent2=exp2, knee=knee_list[k])

        #save peaks to mat file (for later processing with Mat2NlxSpike to generate .ntt file)
        print('saving to mat file')
        outname = OutFolder+"/SimSig_knee"+str(knee_list[k]) +"_"+str(l)+'.mat'
        savemat(outname, {"Timestamps":times, "SimSig": knee_ap, "Params":[exp1,exp2,knee_list[k]]})
        del knee_ap `

[TomDonoghue]

Hey @SusanL82 - I had a check through, and I can't replicate any of the weird spectra from sim_knee.

For example, the following code does the procedure you describe all in Python:

fs = 6000
knee_ap = sim_knee(10, fs, exponent1=-0.5, exponent2=-1, knee=10000)
plot_power_spectra(*compute_spectrum_welch(knee_ap, fs, nperseg=2*fs))

This gets expected power spectra, e.g.:

I explored around different exponent / knee / time / sampling values and didn't see anything suspect. Concatenating segments didn't seem to create an issue, nor did trying some longer sections. I'm not sure what you're seeing - but it may have to do with something outside of the simulation creating, and relate to how you are concatenating / combining / estimating power spectra? If you can share some code that shows this issue directly in Python (for example, adapting the code above), we can revisit - otherwise I can't find anything specific to neurodsp to check / fix here.

Also, sim_knee gets much slower for longer signals, because it computes the expected amplitude per frequency of the FFT. Because the number of FFT frequencies scales with signal length, this will quickly slow down with high sampling rates and/or long signals. It's therefore expected that it takes a long time. We've never really used these approaches for minutes to hours or simulated signals - you may want to consider alternate approaches for simulating such long signals and/or see if you can address your questions with shorter signals.