Discrepancies between Power Spectra and Data

Hi MATLAB users,

I have time-series of repetitive burst, like the following link (lower trace): http://mitpdev.mit.edu/library/erefs/arbib/images/figures/A002_fig001.gif

The bursting occurs at a frequency of 0.5Hz, and within the burst there are frequencies of between 50-70Hz. Sampling at 10millisecs, my power spectra;

Fs=100;
pwelch(event_train,128,120,[],Fs,'onesided')

doesn't give me what I would expect (peaks at 0.5Hz and 50-70Hz). What am I doing wrong?

Thanks,

Linford

[Linford Briant]

So the timestep for my data is 10ms. This may be a confounder here.

Thanks,

L

Answer 1

[Wayne King]

Hi Linford, A sampling frequency of 100 Hz (dt=0.01) is going to result in aliasing for the content from 50-70 Hz. As far as your other question, is there anyway you can post your time series anywhere? How long is your data vector? Your frequency resolution with a segment length of 128 points and your sampling frequency is almost 0.4 Hz so one thing to do would be to increase your segment length.

Without seeing the data or your PSD estimate, it's hard to say anything more. Another thing is what is the mean of the data. Is the mean nonzero?

If so, the power at DC (zero frequency) can obscure the spectral details you are interested in, so you may want to detrend() first.

[Linford Briant]

Hi wayne, thanks for the speedy reply. So I'm slightly confused; I want this second peak for the 50-70Hz peaks to occur. How do I choose a sampling frequency so that it pics up both this low frequency modulation (~0.5Hz) and this higher, broader peak (50-70Hz)? Would I choose one that wasn't at a similar frequency to either of them?

I can get you the datavector, its 420000 in length.

Thanks wayne,

Linford

[Wayne King]

The sampling frequency is determined by the acquisition process, the hardware/software that was used to acquire the spike train. If that process is sampled at 100 Hz, I'm assuming that they use a lowpass filter to prevent aliasing. Those details should be available to you.

For your purposes, You would have to reacquire the data at a higher sampling rate, say 200 Hz to be able to resolve activity in the 50-70 Hz range.

[Linford Briant]

I set the sampling frequency to 0.025ms - the data is from simulations from the neuronal modelling software, NEURON. Obviously with 110second simulations, I need to change the sampling frequency from 40000 to something a bit more managable! So I 'resampled' the data in matlab, to acquire a spike train with 1000Hz sampling (1/ms), and it seems to work fine... I am getting modes of frequencies at 50Hz, and other higher frequencies that I should check exist in my inter-spike-interval histogram. I'm still not getting this 0.5 peak though.

Also, how do I change the y-axis of the output graph from pwelch to something other than power/freq? What would you recommend.

Thanks Wayne,

Linford

[Wayne King]

I see, so this is simulated data? Does the simulation result in data that is zero mean? If you enter mean(event_train) what do you get? pwelch() is a power spectral density estimate so it does give you power/unit frequency. Did you try increasing your segment length? with such a long data vector, you can certainly make that segment pretty long and increase the overlap accordingly.

Can you post the data vector somewhere?

[Linford Briant]

Also, the mean is 0.022... why does a non-0 mean matter? I see some people subtract this from the vector they are computing the power-spectra on. Why is this?

[Linford Briant]

I can email you both the raw-data and the "event-train" that I sample at 200Hz - by converting spikes into delta-functions. What is your email?

[Wayne King]

Search for me on the CSSM newsgroup

http://www.mathworks.com/matlabcentral/newsreader/

[Linford Briant]

Ah-ha, subtracted mean and get a peak at 0.58Hz... perfect!

[Wayne King]

yep, see! That's because when you have a non-zero mean, the value at zero frequency can get very big (you sum all the values in the data). That can obscure other frequency content of interest. :)

Answer 2

[Rick Rosson]

With a sampling rate of 100 samples per second, the Nyquist frequency is 50 hertz. That means the maximum frequency that you can represent in the spectrum is 50 hertz, which is less than your frequencies of interest, 50 to 70 hertz.

[Linford Briant]

Ahhh so I need a sampling frequency of, say 200Hz?

[Rick Rosson]

You need a sampling rate of at least 140 samples per second to see a peak at 70 hertz. So, yes, 200 samples per second would be okay.

[Rick Rosson]

Also, Wayne's suggestions related to the resolution of the spectrum and filtering out the DC component are spot on.

[Linford Briant]

Thanks Rick - most helpful.

Linford

Answer 3

[Dimitris Malamas]

Hi, I have time-series (MER) of a Parkinson patient, the sample frequence is 24000 Hz. I want to estimate the power spectrum density of the spike trains, I used Welch method, it's parameters incuded a Hanning window and a 50% overlap between windows that produced a 0.5 Hz spectral resolution. But i have a problem, i am not sure if i used the right input data, from Matlabs help : [Pxx,f]=pwelch(x,window,noverlap,nfft,fs), i am not sure about fs, it means the sample frequency of MER or something else?