Power Spectra estimation after FFT

Hello guys!

I am working on a frequency analysis for skin conductance. The frequencies expected on those kinds of signals are very low (between 0.045 and 0.25 Hz), so I would like to get the PSD in specific low-frequency ranges, i.e. (0 - 0.045), (0.045 - 0.15), (0.15 - 0.25), and the total power of the signal. To do so, I computed the Fast Fourier Transform, but now the problem is computing the power spectrum after the fourier.

The tutorials from Matlab have not been of much help, and also there are so many functions that can be used that I'm feeling a little lost at the moment! Can anyone help me work out this issue? Among the many options, which one do you think works best, in terms of both efficacy and computational time? I'm not looking for anything too fancy really, just an efficacious and standard protocol.

The processing pipeline I am following is:

1. 8th order Chebyshev Type I lowpass filter, with 0.8 cutoff
2. Downsampling from 4 to 2 Hz (factor = 2)
3. 8th order Butterworth highpass filter with 0.01 cutoff (do you think is appropriate, given that I'm mostly interested in low frequencies?)
4. Welch's periodogram (50% overlap)
5. Blackman window (128 points) applied to each segment
6. FFT computer on each windowed segment

I also attached the code because it would be amazing to have examples or suggestions about how to proceed based on real examples! Any advice is always appreciated!

Thank you very much to any of you who is willing to help!

Luca

[dpb]

W/o data to look at, it's all just a hypothetical exercise with no idea of whether would make sense or not.

The first thing would be to look at what the data collection protocol is for such low frequency signals...takes a long time to get enough data to have decent statistics at 0.05 hz -- that's 20 sec for just one cycle at the high frequency end for your above low range of interest--it only goes up from there.

One would guess just drift and thermal effects, etc. could easily be as much as the actual change in conductance...not to throw a damper on anything, but... :)

[Luca Merolla]

Hello! Thanks for your reply!

I shared my code here, could you please give me some tips? First of all, I hope the filtering and windowing parts are correct. Then my first idea was to compute the area under the curve (AUC) as an estimate of the power spectrum: is this reasonable? When I do the AUC on the FFT of the signal, I get both the real-valued number and its imaginary part: how shuld I interpret this result? Can it still be useful for my purpose?

Are there other methods to extract the PSD from the signal? As I said, I'm not looking for anything too complex, just a standard procedure.

Of course, if you feel like something can be changed and improved please tell me! It's good to have feedbacks :)

Thank you again!

Ps: I know this signal is short (less than 10 minutes), but in the "official" experiment it will be around 30 minutes long

[dpb]

There's a strong linear trend there with a few distinct locations where the trend line is discontinuous (mostly a drop, then continued essentially identical slope). What's the cause for that/

What are the "tag" lines? Is that some stimulus condition? Then there's a huge spike out in the middle not related to anything else noted?

Is your sample rate 4 Hz I gather? Do you have any analog lowpass filtering before the input is sampled? If there is higher energy content in the signal, once it's sampled, it's already there -- post-sampling bandpass filtering won't remove aliasing already extant in the data.

What is the end result for which one is looking? What's the information thought to be in the frequency content of these measurements? (You don't necessarily have to answer that here if it's proprietary in nature, but the basic question of "Why?" is one that always needs answering.)

[Luca Merolla]

No worries, it's no secret!

So this is a skin conductance signal acquired with the Empatica E4 device (sampling frequency is 4 Hz), the tag lines are the beginning and end of a baseline acquisition, and the sharp spikes are most probably artefacts related to electrode drifts. All the filtering is computed offline, as it is a wearable device. The energy content I am expecting to find is very low, in the range of 0.045-0.25 Hz, related to sympathetic activation (like in the low frequency HRV analysis). Indeed, the aim of this analysis is to find features in the signal related to arousal and sympathetic nervous system activity: if that particular frequency band (0.045 - 0.25 Hz) holds the most of the power of the signal, that we can be pretty sure there is a sympathetic activity going on in the subject. This is what I read in the lietrature. Given the features of the signal, a time-domain analysis is not suitable because there isn't much peak activity, so I chose to focus on a frequency analysis.

I hope this can help to understand a little more, if you have more questions I am happy to answer! Thank you very much for your time!

Answer 1

[LO]

try this example:

sampling_freq = 5000
window = 8192, 
noverlap = 4096, 
nfft = 8192,
[p,f] = pwelch(your_signal, 8192,4096,8192,5000);

max_power=10*log10(max(p)); % this would be the max power
plot(f,10*log10(p)); % this would be the power plot

[Luca Merolla]

Hey, thank for the suggestion! Can you please explain why did you choose those values?

[LO]

well they fit my data. is like the binning for the histogram function.

in general nfft is a power of 2 (in my case I use often 2^10), in this case is much higher because I want to estimate the power over a larger sample (but is the sampling freq which determines the time range of what you want to measure the power of... for instance you have a sampling freq of 5000, it means you have 5000 points per sec.

for how many seconds you want to measure of your signal? this is basically the question you need to answer, but there is no rule, it depends on your data and what you want to get out of them.

the window size also, it represents just the amount of data points you want to analyze at each FFt cycle

the overlap improves the "resolution" of such analysis - as it determines the percentage of overlap between adjacent FFT analysis windows (you could set it very low or even very high, normally in my formulas is a fraction of the nfft, such as overlap = round(0.9*nfft) for a 90% overlap).

hope it helps !

[Luca Merolla]

Thanks for the explaination!

So I've been looking at some formulas to calculate an appropriate window length, overlap and nfft, but it seems there is no standard rule. What formula or theory do you usually follow to calculate them? I am interested in the low frequency band of the signal, from 0.045 to 0.25 Hz, and my sampling frequency is 4 (possibly downsampling at 2 Hz will be done, but it's still to be decided if it's appropriate). In general, the pipeline I'm following is:

1. Lowpass filtering (Chebyshev Type I, 8th order, cutoff 0.8 Hz)
2. Downsampling at 2 Hz
3. Highpass filter (Butterworth, 8th order, cutoff 0.01)
4. Calculate power spectra using Welch's periodogram method (50% overlap)
5. Apply Blackman window (128 points) to each segment
6. Calculate FFT for each windowed segment

It would be great if you could give me some indications. Thanks a lot again!

[LO]

Well if my answer was working accept it 😊 For the rest it really depends on your data, the length of the segment you analyze. The frequency of the signal is irrelevant. However if it is a signal with frequency modulations (changes) then you want to use a nfft higher than 2^10 in order to better resolve frequency than time. Conversely if you want a better time resolution go lower. Regarding the other params it depends on the "data points" your signal is composed by. I suggest to get an idea, search for A Wikipedia page or try on your phone a spectrogram app or simply audacity on your computer and see how the different settings affect your spectrogram/fft analysis